AD-A066  155  NORTHROP  RESEARCH  AND  TECHNOLOGY  CENTER  PALOS  VERDES  —ETC  F/G  20/5 

THEORETICAL  MODELING  OF  MOLECULAR  AND  ELECTRON  KINETIC  PROCESSE— ETC <U) 
JAN  79  KB  LACINA  N000I4-78-C-0499 

UNCLASSIFIED  NRTC-79-7R-V0L-1  ML 


£ 


fltTCOPfc  AD  AO  66155 


— — —— — - ll,’l“,ll,“11  ■ 

LEVELS 

THEORETICAL  MODELING  OF  MOLECULAR  AND  ELECTRON 

KINETIC  PROCESSES 

Volume  I 

Theoretical  Formulation  of  Analysis  and 
Description  of  Computer  Programs 

January,  1979 


William  B.  Lacina 


§ 

Sponsored  by 

OFFICE  OF  NAVAL  RESEARCH 
Department  of  the  Navy 

Contract  No.  N000 14-78-C-0499 


NORTHROP 

Research  and  Technology  Center 


79  03  21  016 

j 


THEORETICAL  MODELING  OF  MOLECULAR  ANDJELECTRON 

KINETIC  PROCESSES* 


Volume^ 


JHieoretical  Formulation  of  Analysis  and 

Ascription  of  £pmPuter£ro9ram' 


illiam  B.  Lacina 


NORTHROP  RESEARCH  AND  TECHNOLOGY  CENTER 
One  Research  Park 

Palos  Verdes  Peninsula,  California  90274 
Telephone: (213)  377-4811 


This  documsr 
for  pub'j-  pi 
distribution 


is  unlimited. 


OFFICE  OF  NAVAL  RESEARCH 
Department  of  the  Navy 

Contract  NojN4d4l4-78-C-^499 


1 CCU  At  TV  CU>»»  AlCAnox  0 A TH«*  A ACC  0m 


REPORT  DOCUMENTATION  PAGE 


Lr  rtncnM'!  r f.i 


■LAkMMkMWM 


A.  ACJIACAmihC  2 AC.  *C»Q*T  NUA«CR 

NRTC-79-7R 


CCMTAaCT  2 a cram r hum4 


N000 14-78-C-0499 


/v>^ 


NRTC-79-7R  x 


" THEORETICAL  MODELING  OF  MO-  — Report 

LECULAR  AND  ELECTRON  KINETIC  PROCESSES  '2“  # w|  „ l5  l<?79 
Voi.  I:  Theoretical  Formulation,  of  Analysis  and 
Description  of  Computer  Programs 


William  B.  Lacina  r ' 


*.  ACAAQAANMC  CROAMlUnGN  .*• AMi  amQ  AOCAC33 

NORTHROP  Research  it  Technology  Center  y 
One  Research  Park' 

Palps  Verdes  Peninsula.  Calif.  90274' 


it.  caNraobkiMC  o»nex  ****  mo  «oa»us 

Office  of  Naval  Research,  Dept,  of  the  Navy 
900  N.  Quincy  Street 

Arlington,  Virginia  22217  _ 


14.  MONtrSNIMC  A4SMCY  HAMS  * .eoeuvit  trmm  C« im«M  Q !)>*•) 

Office  of  Naval  Research,  Pasadena,  Calif. 


AMO  A *K3  AK  MUM4CAS 

Attn. : 613C:MAK 


•a.  oats 

January,  1979 


UNCLASSIFIED 


AtiUTtOf#  ST  JkT"Vd€)iT  (mi  tut  a Aim  met) 

I DISTRIBUTION  STATEMENT  A 


Unlimited. 


Approved  hot  public  xeleasei 
Distribution  Unlimited 


it.  strnieurtQM  r »*  .<«•«  m *>««»  :o.  ■/  tram 


nnTTymi 


tf.  ACT  (Csinww  impm«  *««»  </  M t4mmnty  *y  »<««* 

Molecular  kinetics,  electron  kinetics,  Boltzmann  equation,  laser  kinetics, 
transient  analysis,  electrically  excited  laser  KrF  laser,  rare  gas  halides 
kinetic  modeling,  superelastic  collisions,  electron- electron  collisions 


A CAT  A ACT  (Camummm  am  ffWf  tlO  »f  liantity  ly  tl««A 

This  report  describes  and  documents  a comprehensive  and  reasonably  general 
computer  analysis  applicable  to  a broad  class  of  transient,  electrically  excited 
laser  systems.  The  theoretical  model  is  formulated,  for  a spatially  homogene- 
ous medium,  in  terms  of  the  coupled  set  of  equations  which  describe  the  mole- 
cular kinetics,  electron  kinetics,  external  discharge  circuit,  and  optical  radia- 
tive extraction.  One  of  the  unique  features  of  the  present  code  is  that  the  mole- 
cular kinetic  subroutines  for  an  arbitrary  reaction  scheme  are  synthesized  auto 
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matically  by  translation  of  an  input  queue  of  symbolically  expressed  reactions 
into  computer- coded  FORTRAN  equations,  which  are  automatically  coupled  to 
a Boltzmann  analysis  of  electron  kinetics.^ 

The  molecular  kinetics  reaction  scheme  is  translated  into  subroutines  which  pro- 
duce the  instantaneous  rates  n*  and  the  Jacobian  Sr^/Sn.  (for  the  population  den- 
sities and  electrical  circuit  parameters)  so  that  a multistep  Gear  technique  can 
be  implemented  for  integration  of  the  ("stiff')  system  of  coupled  equations.  The 
radiation  fields  for  laser  extraction  are  formulated  in  terms  of  the  intracavity 
photon  density,  with  a source  term  from  spontaneous  noise  and  amplification  (and 
absorption)  by  the  excited  medium.  This  approach  allows  radiation,  electrons, 
and  all  of  the  other  molecular  "species"  to  be  treated  on  a parallel  basis. 

-^The  electron  kinetics  analysis  consist's  of  numerical  solution  of  the  Boltzmann 
equation  for  the  electron  energy  distribution.  The  technique  for  solution  per- 
mits inclusion  of  inelastic  binary  electron- molecule  collisions,  superelastic 
collisions,  electron- electron  collisions,  momentum  transfer  (with  recoil),  a 
source  term  .Sq-S-fw)  1*  6c  s(u^for  external  creation  of  electrons  (at  zero  energy 
and  over  a distribution  of  energies)  by  sources  of  ionization,  and  retention  of 
the  term  proportional  to  dne/dt  for  situations  involving  a net  change  in  the  elec- 
tron density.  From  the  electron  energy  distribution,  all  of  the  plasma  parame- 
ters, electron  excitation  rates,  and  power  partitioning  can.be  obtained.  ^ 

The  general  laser  kinetics  program  described  has  been  constructed  with  consi- 
derable emphasis  on  flexibility  and  simplicity  of  usage  for  the  user.  In  addition 
to  the  instructions  for  usage  described  in  the  present  report,  the  FORTRAN 
source  decks  contain  extensive  internal  COMMENT  card  documentation  as  well. 
Even  synthesized  subroutines  which  are  generated  by  translation  of  the  input  re- 
action queue  contain  COMMENT  card  documentation.  In  addition  to  the  laser 
kinetics  synthesis  and  analysis  code,  an  independent  program  for  solution  of  the 
Boltzmann  equation  has  been  developed,  and  is  also  included  in  the  present  re- 
port. All  of  the  programs  and  subroutines  have  been  written  in  FORTRAN  IV  for 
Extended  Fortran  Compilation  on  the  CD<2  6000  and  CYBER  Series  of  computers 
(which  use  a 60- bit,  10- BCD  display  character)  word.  Extensive  modifications 
may  be  required  if  these  codes  are  to  be  executed  on  some  other  system. 

The  present  programs  were  developed  in  support  of  experimental  research  pro- 
jects for  the  KrF  excimer  laser.  Therefore,  input  and  output  for  analysis  of  a 
KrF  system  are  used  for  illustration.  Program  listings  are  given  in  Vol.  II  of 
this  report. 
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1.0  INTRODUCTION 


Northrop  Research  and  Technology  Center  has  been  involved  in  a variety 
of  research  programs  which  have  resulted  in  significant  contributions  to 
the  development  of  high  power  laser  technology  as  well  as  advances  in  the 
understanding  of  fundamental  laser  physics.  Recent  investigations  have 
been  focused  on  uv  and  visible  laser  concepts,  with  primary  emphasis  on 
the  rare  gas  halide  excimer  systems  such  as  KrF  . In  theoretical  sup- 
port of  these  experimental  investigations,  the  present  author  has  devel- 
oped a comprehensive  computer  code  which  is  sufficiently  general  that  it 
would  be  of  substantial  future  benefit  for  a wide  variety  of  analytical  stud- 
ies of  lasers,  chemical  kinetics,  and  electric  discharge  phenomena. 

This  code,  which  contains  far  more  generality  and  flexibility  than  any  of 
its  nominal  competitors  currently  in  existence,  will  be  thoroughly  de- 
scribed and  documented  in  the  present  report  in  order  to  make  it  avail- 
able to  the  general  community. 

During  the  past  several  years,  there  has  been  a continuing  interest  and 
extensive  effort  in  theoretical  modeling  of  laser  kinetic  phenomena.  In 
support  of  research  directed  toward  development  of  new  laser  concepts, 
construction  of  a computer  analysis  to  model  the  fundamental  physical 
mechanisms  and  microscopic  kinetic  processes  can  provide  a valuable 
tool  for  understanding  the  laser  operation  and  for  optimizing  its  perfor- 
mance. Such  a model  is  continuously  refined  as  understanding  of  the  la- 
ser reaction  scheme  evolves. 

Exploratory  research  for  the  development  of  promising  new  laser  systems 
is  often  a combination  of  art  and  science,  and  typically  requires  both  ex- 


* This  code  was  developed  under  internal  IR&D  funding,  with  partial 
support  from  the  Advanced  Research  Projects  Agency  of  the  Depart- 
ment of  Defense,  monitored  by  the  Office  of  Naval  Research  under 
contract  N00014-76-C- 1100. 
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perimental  and  theoretical  understanding  of  several  overlapping  disci- 
plines such  as  molecular  spectroscopy  and  kinetics,  electron  kinetics 
and  discharge  circuitry,  fluid  dynamics,  and  laser  radiative  extraction. 
The  first  step  in  the  discovery  of  a new  laser  system  is  the  formulation 
of  a conceptual  mechanism  of  excitation  and  energy  transfer  for  a set  of 
candidate  atomic  and  molecular  species;  usually,  this  is  based  initially 
upon  only  am  intuitive  or  semiquantitative  understanding  of  the  spectro- 
scopy and  kinetic  processes  for  the  proposed  system.  Inspiration  for  a 
concept  often  originates  from  current  fundamental,  research  in  the  spec- 
troscopy of  highly  excited  states  of  known  molecules  or  unfamiliar  trans- 
ient species,  and  from  studies  of  collision  dynamics  and  measurements 
of  energy  transfer  reaction  rates.  If  initial  experimental  investigations 
confirm  that  a hypothetical  laser  concept  can  be  successfully  realized, 
more  detailed  studies  complemented  by  theoretical  modeling  and  analysis 
are  required  in  order  to  optimize  laser  performance  and  to  formulate 
quantitative  scaling  laws  for  extrapolation  to  high  power,  efficiency,  and 
volume.  For  this  purpose,  it  is  necessary  to  understand  the  kinetic  and 
pumping  processes  which  are  responsible  for  the  formation  and  quenching 
of  the  excited  species.  A theoretical  model  of  the  fundamental  physical 
mechanisms,  implemented  by  a computer  code,  is  a useful  asset  for  the 
development  and  optimization  of  a practical  laser  device. 

The  typical  development  of  a laser  kinetic  model  proceeds  by  construct- 
ing an  analysis  of  a spatially  uniform  gain  medium  in  a simple  plane  par- 
allel resonator  with  an  attempt  to  describe,  as  completely  as  possible, 
all  of  the  coupled  phenomena  involving  atomic  and  molecular  kinetics, 
electron  kinetics,  discharge  circuitry,  optical  radiative  extraction,  fluid 
dynamics,  etc.  which  characterize  the  laser  system.  Although  the  goal 
of  such  an  effort  is  to  provide  a reliable  analytical  tool  for  optimizing 
the  development  of  a practical  laser  device,  such  a model  may  also  be 
useful  even  when  the  laser  may  fail  to  accurately  exhibit  the  predicted 
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behavior  in  all  respects.  In  such  a case,  one  may  still  be  able  to  define 
parametric  trends  which  suggest  the  necessary  approach  for  optimization 
of  the  laser  device.  Even  with  an  incomplete  model,  it  is  often  possible 
to  better  understand  the  importance  of  microscopic  processes  or  excited 
species,  and/or  to  eliminate  some  unimportant  factors  from  further  in- 
vestigation. In  fact,  evolution  in  the  theoretical  understanding  of  the  fun- 
damental physical  mechanisms  of  the  laser  proceeds  by  just  such  a pro- 
cess in  which  the  sensitivity  of  various  r^te  constants  is  assessed,  the 
dependence  of  observed  output  on  various  experimental  parameters  is  de- 
termined, and  the  kinetic  reaction  scheme  is  modified  to  account  for  and 
to  better  explain  the  observed  phenomena.  Although  such  a model  is  con- 
tinually refined  as  the  understanding  of  the  reaction  scheme  evolves,  even 
preliminary  predictions  from  the  analysis  may  be  sufficiently  reliable  to 
justify  more  extensive  experimental  effort,  or  to  define  and  narrow  the 
range  of  parameters  or  direction  of  approach. 

The  most  important  theoretical  components  of  the  model  consist  of  the 
analysis  of  the  fundamental  physical  mechanisms,  and  in  particular,  the 
refinement  of  the  molecular  kinetic  reaction  scheme  and  estimates  of  its 
rate  constants.  Simultaneously,  attempts  at  more  sophisticated  descrip- 
tion of  related  problems  such  as  discharge  stability,  medium  nonuniform- 
ity, optical  extraction  and  resonator  configuration,  flow  properties,  and 
mode-medium  interactions  can  be  undertaken.  In  general,  all  of  the  lat- 
ter refinements  are  concerned  more  with  an  improved  description  of  ef- 
fects which  result  from  spatial  inhomogeneity  rather  than  with  fundamen- 
tal physical  mechanisms,  and  are  therefore  most  appropriate  only  in  the 
final  stages  of  development  of  a realistic  engineering  model.  That  is 
not  to  say  that  these  effects  are  unimportant,  for  they  often  present  prac- 
tical limitations  that  make  it  difficult  to  attain  the  results  predicted  for 
an  idealized  homogeneous  medium. 
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Therefore,  the  scope  of  the  present  work  was  to  develop  a coupled  analy 
sis  of  molecular  kinetics,  electron  kinetics,  electric  discharge  circuit, 
and  optical  radiative  extraction,  which  form  the  fundamental  basis  for 
the  description  of  the  physics  of  an  electrically  excited  laser  system. 
The  goal  of  any  such  model  is  to  successfully  predict,  for  a spatially 
uniform  gain  medium,  the  excited  state  population  densities,  gain  and 
saturation  parameters,  power  transfer  and  extraction  rates,  electrical 
power  partitioning,  optical  conversion  efficiency  and  specific  power,  op- 
tical output  power  and  spectral  distribution,  etc. 


in  general,  the  type  of  model  developed,  the  detail  and  approximations 
considered,  the  computational  sophistication,  the  input/output  flexibility 
of  the  computer  code,  etc.,  can  vary  considerably  depending  upon  the  spe 
cific  problem  and  application.  There  is  currently  much  needless  dupli- 
cation of  effort  devoted  to  laser  kinetic  modeling.  For  any  specific  laser 
concept,  there  is  often  a considerable  overlap  as  well  as  a frequent  lack 
of  communication  between  competing  investigators  engaged  in  identical 
kinetic  analyses.  Furthermore,  this  overlap  is  even  more  extensive 
when  one  considers  the  work  of  those  who  have  developed  essentially  the 
same  analytical  models  for  different  systems  which  (although  not  identi- 
cal) are  nevertheless  members  of  a broad  class  of  electrically  excited 
lasers  for  which  a more  cenerallv  andicable  analysis  mni ri  h.va 
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al  times,  and  future  programs  cannot  directly  build  upon  or  benefit  from 
the  previous  work.  It  is  impractical  to  construct  a new  computer  code 
for  the  analysis  of  each  new  kinetic  reaction  scheme.  If  the  necessity  of 
computer  programming  were  substantially  eliminated  by  the  availability 
of  an  applicable,  user-oriented,  existing  computer  code,  a significant 
amount  of  time  could  be  saved  and  allocated  to  more  productive  efforts 
to  understand  the  physics  of  the  laser  system. 

Recognition  of  the  lack  of  a coherent  and  unified  approach  to  the  problem 
of  laser  kinetic  modeling  provided  the  inspiration  to  undertake  the  devel- 
opment of  a generalized  computer  code  which  would  be  applicable  to  a 
broad  class  of  transient,  electrically  excited  laser  systems.  The  pre- 
sent code  automatically  synthesizes,  for  an  arbitrary  reaction  scheme, 
a completely  self  contained  computer  code  for  numerical  solution  of  the 
coupled  equations  which  describe  the  molecular  kinetics,  electron  kine- 
tics, external  driving  circuit,  and  optical  radiative  extraction.  This 
code  automatically  generates  its  own  subroutines  for  the  analysis  of  mo- 
lecular kinetics  by  translating  symbolic  reactions,  provided  as  input, 
into  computer- coded  equations.  Thus,  for  the  most  complicated  reaction 
scheme  containing  an  arbitrary  number  of  kinetic  collision  processes  and 
interacting  species,  it  is  possible  to  obtain  the  complete  computer  code 
required  with  virtually  no  effort.  The  automatically  synthesized  molecu- 
kinetics  subroutines  are  combined  with  a master  executive  program 
and  all  other  required  subroutines  (e.g.  for  numerical  solution  of  the 
Boltzmann  transport  equation)  to  form  a completely  self-contained  cou- 
pled analysis  based  upon  the  specified  reaction  scheme. 

With  appropriate  modifications,  this  code  can  be  easily  adapted  to  the 
analysis  of  related  problems  in  various  fields  which  involve  studies  of 
chemical  kinetics.  Its  internal  subroutines  for  numerical  solution  of 
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the  Boltzmann  equation  can  be  employed  for  the  analysis  of  transport 
properties  of  electrical  discharges,  and  a separate  program  developed 
for  this  purpose  will  be  described  in  Sec.  3.4  (cf.  Vol.  II  for  listing). 

The  present  computer  code  is  completely  user-oriented,  and  for  a very 
broad  class  of  transient,  electrically  excited  laser  systems,  it  complete- 
ly eliminates  the  necessity  for  the  laborious  task  of  constructing  new 
computer  codes.  It  will  do,  automatically,  work  that  would  require  sev- 
eral months  if  done  manually,  and  it  will  provide  an  input/output  flexibi- 
lity for  the  synthesized  analysis  surpassing  that  which  typically  charac- 
terizes existing  codes  that  have  been  developed  for  specific  systems. 

It  was  felt  that  documentation  and  dissemination  of  this  computer  code  to 
those  engaged  in  laser  kinetic  modeling  would  provide  substantial  bene- 
fit and  economy  by  the  elimination  of  needless  repetition  and  duplication 
of  effort.  It  has  been  the  intent  of  the  author  that  every  effort  be  made 
to  maximize  the  flexibility  and  simplicity  of  usage  of  the  code,  and  that 
the  description  of  the  analysis  and  computational  structure  be  sufficiently 
clear  and  detailed  that  others  shall  be  able  to  modify  and/or  extend  the 
analysis  as  they  require.  The  objective  has  been  to  reduce  mechanical 
programming  considerations  to  a minimum,  and  to  provide  the  user  with 
(limited)  diagnostic  assistance  as  well  as  some  protection  by  automatic 
exit  if  specified  error  conditions  are  encountered  during  execution. 

The  present  code  has  been  written  in  an  extended  version  of  FORTRAN 
IV  with  syntax  and  conventions  compatible  with  the  CDC  6600  computer. 

It  makes  extensive  use  of  a 60-bit,  10- BCD  character  word  size,  as  well 
as  various  special  features  permited  by  the  current  CDC  FTN  compiler, 
so  considerable  effort  may  be  required  if  it  should  be  necessary  to  adapt 
it  to  execute  on  other  machines.  However,  because  the  CDC  6600  cur- 
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rently  enjoys  widespread  usage  for  scientific  applications  at  government, 
industrial,  and  university  research  laboratories,  it  is  anticipated  that 
the  code  will  be  suitable  for  the  majority  of  potential  users.  In  order 
that  the  program  may  be  more  easily  modified  or  extended,  the  source 
decks  have  been  extensively  documented  with  internal  COMMENT  cards. 
Furthermore,  even  the  synthesized  subroutines  are  generated  with  inter- 
nal COMMENT  card  documentation  to  make  them  completely  readable  and 
understandable.  The  synthesized  code  is  a complete  and  self-contained 
FORTRAN  source  program. 

Obviously,  because  this  code  has  been  developed  as  a research  tool  in 
support  of  current  experimental  laser  development  programs,  it  is  not 
possible  to  guarantee  that  it  will  perform  perfectly  under  all  conditions 
to  which  it  may  be  applied,  although  every  effort  has  been  taken  to  in- 
sure that  it  do  so.  Even  programs  which  are  free  of  logical  or  program- 
ming  errors  can  fail  for  reasons  of  numerical  ill- conditioning,  converg- 
ence difficulties,  computational  instabilities,  accidental  division  by  zero, 
etc.  Because  the  present  code  was  developed  primarily  in  support  of 
KrF  laser  research,  the  illustrative  examples  to  be  presented  in  subse- 
quent sections  shall  apply  to  that  system. 

It  has  been  said  that  "All  computer  programs  must  be  assumed  to  con- 
tain bugs  until  proven  otherwise,  which  is  impossible."  Users  who  dis- 

1 cover  suspected  "bugs"  in  the  present  analysis  are  encouraged  to  com- 

municate them  to  the  author  who,  in  turn,  will  attempt  to  offer  assistance 
to  users  when  necessary.  It  is  hoped  that  investigators  who  benefit  from 
the  availability  of  this  code  will  make  appropriate  acknowledgments  of 
its  usage  in  the  footnotes  or  references  of  their  published  papers,  con- 
ference presentations,  or  reports. 
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The  scope  of  the  present  work  was  the  development  of  a comprehensive 
and  reasonably  general  computer  analysis  applicable  to  a broad  class  of 
transient,  electrically  excited  laser  systems.  The  theoretical  model  is 
formulated,  for  a spatially  homogeneous  medium,  in  terms  of  the  coupled 
set  of  equations  which  describe  the  molecular  kinetics,  electron  kinetics, 
external  discharge  circuit,  and  optical  radiative  extraction.  For  each 
new  laser  kinetic  scheme,  the  equations  which  define  the  molecular  and 
electron  kinetics  are  different,  and  the  main  advantage  of  the  present 
computer  code  is  that  it  generates  the  required  analysis  automatically  by 
synthesizing  the  necessary  subroutines.  It  is  especially  useful  for  sys- 
tems characterized  by  complicated,  but  relatively  unsystematic,  kinetic 
reaction  schemes.  Its  most  appropriate  application  is  to  uv  and  visible 
laser  systems  which  involve  excitation,  energy  transfer,  and  stimulated 
emission  from  electronic  states  of  atoms  and  molecules. 

For  infrared  laser  such  as  CO,  CO^,  HF,  etc.,  the  mechanism  of  pump- 
inS»  energy  transfer,  and  radiative  emission  involves  the  excited  vibra- 
tional levels  of  the  ground  electronic  state,  and  for  these  lasers,  the  pre- 
sent approach  (with  modifications)  would  be  relatively  inefficient,  if  not 
totally  inadequate.  For  the  infrared  laser  systems,  a more  complicated 
treatment  of  radiation  must  often  be  included  to  allow  for  possible  absorp- 
tion and  emission  for  several  vibrational- rotational  transitions,  and  there 
are  strong  temperature  effects  in  the  vibrational  kinetics  as  well  as  in  the 
radiative  processes.  In  its  present  form,  the  analysis  does  not  include 
effects  of  changes  in  the  molecular  temperature,  and  there  is  no  provi- 
sion for  temperature  dependent  parameters.  For  such  problems,  the 
construction  of  a more  specialized  analysis  which  explicitly  recognizes 
the  systematic  structure  of  the  vibrational- rotational  level  interactions, 


8 


and  which  can  most  easily  incorporate  temperature  dependent  effects  or 
other  unique  features,  would  probably  be  more  suitable  than  attempts  to 
modify,  extend,  and  apply  the  present  program.  It  is,  of  course,  antici- 
pated that  there  will  be  situations  for  which  extensions  or  modifications 
of  the  present  code,  or  of  the  subroutines  which  it  synthesizes,  will  be 
the  optimum  approach. 

The  components  of  the  present  model,  labeled  with  some  brief  descrip- 
tive remarks,  are  summarized  in  Fig.  2. 1.  From  the  more  detailed 
discussion  of  the  physics  to  be  given  below,  it  should  become  more  ap- 
parent for  which  problems  the  present  analysis  is  most  applicable.  A 
discussion  of  the  limitations,  approximations,  numerical  techniques,  and 
range  of  applicability  of  the  analysis  will  be  provided. 

2.  1 Molecular  Kinetics 

The  molecular  kinetics  are  described  by  the  master  equation  for  the  pop- 
ulation densities  of  all  the  atomic  and  molecular  species  present  in  the 
electrically  excited  gas  mixture,  including  electrons,  ions,  neutrals,  ex- 
cited states,  excimers,  etc.  For  a given  set  of  reactions,  the  formal 
construction  of  the  master  equation  for  the  population  densities  is  illus- 
trated in  Fig.  2.  2.  If  the  species  involved  are  X , X , X , ....  X , any 
reaction  in  the  molecular  kinetic  scheme  can  be  written 

Vl  + "2X2  + »*JX3  + •••  «nX„  Z=£L 

kr 

UlXl  + VZX2  * UiXi  * ■■■  *vn\  +AE  <» 

where  kf  and  k_  are  forward  and  reverse  rate  constants  (s'1,  cmJ/l, 

()  . 

cm  Is,  etc.),  and  and  V.  are  {nonnegative)  integer  coefficients  (most 
of  which  are  zero  or  one).  The  conservation  of  energy  in  Eq.  (1)  gives 
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• MOLECULAR  KINETICS 

Master  Equation  for  population  densities  of  all  atomic  and  molecular 
species  (electrons,  ions,  neutrals,  excited  states,  excimers,  etc. ). 

• ELECTRON  KINETICS 

Boltzmann  transport  equation  for  the  electron  energy  distribution 
to  determine  plasma  parameters,  electrical  excitation  rates,  power 
partitioning,  and  discharge  efficiency. 

• CIRCUIT  ANALYSIS 

External  circuit  equations  to  describe  the  electrical  power  loading  in 
a low  impedance  gas  discharge,  whose  conductivity  is  function  of  L 

• OPTICAL  EXTRACTION 

Oscillator  analysis,  formulated  in  terms  of  intracavity  photon  density, 
with  transient  build-up  of  laser  mode  from  spontaneous  emission  and 
subsequent  relaxation  to  quasi-steady  condition  with  gain  at  threshold 


Fig.  2.1:  Basic  components  for  the  theoretical  analysis  of  an 
electrically  excited  laser  system. 
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MOLECULAR  KINETICS 


Fig.  2.3:  Construction  of  the  molecular  kinetics  master  equation  for  the 
rates  of  change  of  population  densities  for  an  arbitrary  number 
of  interacting  species  in  a general  reaction  scheme  with  colli- 
sion processes  a. 


an  expression  for  the  net  energy  change  AE: 


r 

I 


S MiEi  = 2 yiEi  + AE 

i i 

-3  - 1 

The  net  rate  per  unit  volume  (cm  s ) R for  the  forward  and  reverse 
reaction  of  Eq.  (1)  is  given  by 

R = k£[Xl]"l[X2]“2[X3],‘3  ,..[X/” 

- kr  [Xx]  2[X3]  3...[Xn]  (2) 

where  [X^J  is  the  population  density  (cm  of  the  ith  species  X^.  The 
contribution  that  such  a reaction  makes  to  the  rate  of  production  (or  loss) 
of  the  ith  species  is  given  by 

d/dt  [X.]  = ( Vi  - M.)R.  (3) 

Thus,  for  an  arbitrary  reaction  scheme  containing  several  collision  pro- 
cesses (labeled  by  a),  the  complete  master  equation  becomes 


d/dt  [X.]  = £(>.(«)-  M.(a)]R(cr)  (4) 

a 

If  Eq.  (4)  is  multiplied  by  the  energy  E^  of  the  ith  species  X^,  and  then 
summed  over  i,  an  equation  for  the  conservation  of  energy  (in  terms  of 
power  balance)  is  obtained: 
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d/dt  ^ E.[X.]  = - ^R(a)AE)a) 


(5) 


where 


AE(a)  = ^ [/*.(<*)  - w.(a)]  E. 


(6) 


is  the  net  energy  change  for  reaction  a.  The  sum  over  i in  Eq.  (5), 

(6)  implicitly  excludes  secondary  electrons  e (u),  whose  contribution 
to  the  total  power  conservation  equation  must  be  expressed  as  an  inte- 
gration over  the  continuum  of  electron  energies  u.  The  formulation  of 
electrical  power  balance  will  be  given  in  Sec.  2.  3 below,  where  analy- 
sis of  electron  kinetics  based  upon  the  Boltzmann  equation  for  the  elec- 
tron energy  distribution  is  discussed.  The  total  power  balance  equa- 


tion will  be  developed  in  Sec.  2.  5,  where  the  physical  interpretation  of 
various  terms  corresponding  to  optical,  electrical,  and  kinetic  heating 
mechanisms  will  be  given. 


We  may,  for  convenience,  formally  define  the  parameter  E^  to  be  zero 
for  secondary  electrons,  so  that  their  exclusion  from  the  sums  in  Eq. 
(5)  and  (6)  does  not  need  to  be  made  explicit.  To  be  consistent,  it  is 
then  necessary  (for  the  purpose  of  Eq.  (5))  to  suppress  the  energy  de- 
pendence denoting  the  continuum  of  electron  species  e (u).  In  doing 
this,  it  should  be  emphasized  that  the  parameter  E^  = 0 for  the  single 
species  "electron"  has  no  physical  significance,  and  that  the  power  bal- 
ance relation  (5)  must  ultimately  be  supplemented  with  the  contribu- 
tions originating  from  the  electrical  processes. 


To  write  the  subroutines  which  define  the  molecular  kinetics  for  a com- 
plicated reaction  scheme  is,  in  general,  a difficult  and  time-consuming 
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task,  and  the  resulting  computer  program  would  have  little  flexibility 
for  analysis  of  any  system  except  those  in  a very  limited  class.  (For 
example,  there  have  been  as  many  as  ~ 80  reactions  considered  to  be 
of  possible  importance  for  modeling  KrF  laser  kinetics.)  Furthermore, 
as  understanding  of  a laser  reaction  scheme  evolves,  addition  of  new 
reactions  (or  deletion  of  old  ones)  may  be  continually  required,  in  addi- 
tion to  making  simple  updated  estimates  of  the  rate  constants  for  the 
reactions  retained.  Thus,  such  a code  would  itself  have  to  be  modified 
in  a continuing  manner,  rather  than  merely  executed  with  revised  input 
values  for  the  rate  constants. 

Therefore,  in  order  to  provide  a more  powerful  analytical  capability 
applicable  to  a wide  class  of  problems,  a generalized  code  has  been  de- 
veloped which  automatically  synthesizes  the  coupled  analysis  described 
in  Fig.  2.  1 for  an  arbitrary  reaction  scheme.  The  input  to  this  code  is 
a sequence  of  reactions  and  (initial)  estimates  of  rate  constants.  The 
syntax  for  the  reactions  is  very  flexible,  with  a free  format  that  accepts 
the  reaction  just  as  it  would  normally  be  written.  The  content  of  each 
reaction  is  analyzed,  and  the  appearance  of  each  new  species  is  recog- 
nized. The  syntax  of  each  reaction  is  subjected  to  numerous  tests  to  de- 
tect errors.  If  the  reaction  is  determined  to  be  acceptable,  it  is  decom- 
posed and  translated  into  appropriate  computer- coded  expressions  in 
order  to  generate  synthetically  the  required  subroutines.  Program  gen- 
eration and  execution  are  protected  by  diagnostic  assistance  and  auto- 
matic exit  (according  to  input  requests)  if  specified  error  conditions  are 
encountered.  A more  detailed  discussion  of  the  structure  of  the  code 
will  be  given  in  a subsequent  section  devoted  specifically  to  program  des- 
cription. 

Fig.  2.3  presents  a schematic  flow  diagram  of  the  present  approach. 

The  rate  constants  initially  assumed  in  the  generation  of  the  program 
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can  be  changed  in  the  subsequent  execution,  if  desired.  The  main  pur- 
pose of  the  rate  constants  in  the  initial  input  deck  is  to  define  whether 
the  forward  and/or  reverse  process  is  to  be  included  in  the  rate  expres- 
sions. If,  during  program  generation,  a zero  rate  constant  is  encoun- 
tered for  any  process  (other  than  secondary  electron  collisions),  no 
translation  of  the  forward  (or  backward)  term  occurs  in  the  synthesized 
subroutines.  Of  course,  if  the  correct  values  of  the  rate  constants  are 
known,  they  may  be  entered  for  once  and  for  all  in  the  original  input 
deck.  As  Fig.  2.  3 indicates,  the  synthesized  molecular  kinetics  sub- 
routines are  automatically  combined  with  the  master  executive  program 
and  all  other  required  subroutines  to  form  a completely  self-contained 
coupled  analysis  based  upon  the  specified  reaction  scheme.  The  section 
enclosed  in  the  upper  box  in  Fig.  2.3  is  required  only  when  a new  reac- 
tion scheme  is  introduced  and  generation  of  a new  program  is  required. 
For  all  subsequent  calculations  based  upon  the  given  reaction  scheme, 
the  initial  conditions,  experimental  parameters,  rate  modifications  (if 
any),  and  I/O  and  control  parameters  are  entered,  and  the  complete 
code  constructed  in  the  upper  box  is  executed.  That  is,  only  the  lower 
section  of  Fig.  2.3  is  involved  thereafter. 

There  are  several  obvious  advantages  to  this  approach.  First  of  all, 
there  is  the  simplicity  of  automation  and  the  minimization  of  the  possi- 
bility of  error  in  the  writing  of  complicated  subroutines.  Secondly,  the 
program  diagnoses  error  conditions  in  the  reaction  syntax  which  may 
not  have  been  noticeable  or  which  may  have  been  overlooked.  For  exam- 
ple, duplicate  reactions  are  detected  (even  when  written  backwards), 
charge  particle  and  heavy  particle  conservation  is  insured,  detail  bal- 
ance relations  for  binary  collision  processes  are  enforced,  and  miscel- 
laneous other  error  conditions  are  detected.  Secondary  electron  colli- 
sion processes  are  recognized  and  properly  coupled  to  the  electron  kine- 


GENERALIZED  SYNTHESIS  LASER  KINETICS  CODE 


Synthesis : 


Fig.  2.3:  Schematic  flow  diagram  of  the  functional  structure  of  the 
laser  kinetics  code.  A complete  FORTRAN  source  pro- 
gram is  synthesized  by  translation  of  symbolic  reactions 
into  computer- coded  equations.  Execution  of  the  analysis 
occurs  in  the  lower  box. 
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tics  analysis,  except  for  cases  without  electric  discharge  (e.g.  e -beam 
excitation  only),  in  which  case  their  rate  constants  can  be  defined  by  in- 
put. Finally,  the  subroutines  are  constructed  in  such  a way  that  null  op- 
erations (multiplications  or  additions  involving  zero)  are  completely  eli- 
minated (if  no  rate  constant  is  provided  for  the  forward  or  reverse  pro- 
cess, no  translation  of  the  corresponding  term  occurs),  and  in  such  a 
way  as  to  optimize  computational  efficiency  (repetitious  or  unnecessary 
multiplications  are  minimized).  Even  for  relatively  simple  problems 
with  uncomplicated  molecular  kinetic  reaction  schemes,  the  availability 
of  a program  which  can  automatically  provide  the  investigator  with  a 
user- oriented  code  with  good  I/O  flexibility  and  simplicity  of  usage  is 
of  obvious  advantage. 

It  is  usually  the  case  that  the  rate  constants  involved  in  molecular  kine- 
tic reaction  schemes  can  vary  over  several  orders  of  magnitude,  and 
therefore,  the  resulting  master  equation  (4)  often  becomes  a "stiff"  sys- 
tem of  differential  equations.  Therefore,  the  approach  which  has  been 
taken  for  integration  of  these  equations  is  to  employ  a multistep  tech- 
nique developed  by  Gear^.  This  method  automatically  adjusts  the  inte- 
gration step  size  as  the  solution  proceeds,  in  such  a way  that  required 
accuracy  conditions  are  maintained.  The  Gear  method  requires  subrou- 
tines not  only  for  the  rates  of  change  m = d [X^  /dt  of  the  population 
densities,  but  also  for  the  Jacobian,  6m  / 6m,  as  a function  of  time. 
Because  numerical  evaluation  of  the  Jacobian  is  not  generally  satisfac- 
tory, it  is  necessary  to  generate  both  such  subroutines  symbolically  in 
the  synthesis  section  of  the  program,  where  the  reaction  scheme  is  de- 
composed and  translated  into  computer- coded  equations. 

The  entry  of  control  parameters,  experimental  parameters,  revised  rate 
constants,  initial  conditions,  output  requests,  etc.,  is  quite  flexible,  and 


permits  the  code  to  be  executed  for  a variety  of  situations  of  interest. 

For  example,  rate  constants  for  secondary  electron  collisions  are  ob- 
tained self- consistently  from  the  coupled  Boltzmann  analysis  for  the  case 
of  an  electric  discharge.  For  the  case  of  e -beam  excitation  only,  these 
rate  constants  default  automatically  to  zero,  but  can  be  specified  by  in- 
put, if  desired.  This  is  useful  for  entering  thermal  (i.  e.  room  tempera- 
ture) values  for  electron  recombination  or  attachment  processes,  for 
example.  Thus,  the  same  general  code  can  be  used  for  both  discharge 
or  e -beam  excitation  conditions.  Likewise,  it  is  possible  to  apply  it 
to  either  an  oscillator  or  amplifier  analysis. 

The  integration  of  the  coupled  set  of  equations  over  the  total  specified 
pulse  length  is  carried  out  with  the  Gear  technique,  which  automatically 
adjusts  the  step  size  to  maintain  accuracy  and  stability.  However,  the 
total  pulse  length  is  divided  into  a finite  number  (nominally  50)  of  subin- 
tervals, at  which  time  the  electron  kinetics  are  updated  and  a variety  of 
output  option  requests  can  be  supplied. 


2.  2 Radiative  Extraction 

The  basic  features  of  the  radiation  analysis  are  summarized  in  Fig.  2.4. 
Because  considerable  simplification  can  be  achieved  by  neglecting  com- 
plications associated  with  a rigorous  description  of  the  optical  fields  in 
a realistic  resonator  configuration,  a spatially  homogeneous  medium  is 
assumed,  with  an  optical  resonator  formed  by  two  plane  parallel  mirrors. 
The  net  amplification  of  the  intracavity  radiation  intensity  I(U)  from  the 
stimulated  emission  and  absorption  process 

X*  + 1(1/)  -y— X + I(V),  (7) 
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OPTICAL  EXTRACTION 


• Spatially  uniform  intracavity  optical  fields. 

• Plane  parallel  optical  resonator  (100%,  R),  separation  l_c. 

• Spatially  uniform  gain  medium  of  length  Lg. 

• Distributed  output  coupling  and  loss  (Y)  coefficients. 

• Cavity  threshhold  gain  coefficient: 

ath  • [Y  + ( 1/2)  In (1  / R ) 1 / L_ 

• Buildup  of  optical  fields: 

Spontaneous  emission  in  solid  angle  Q defined  by  output 
coupling  aperture  and  resonator  length;  A/L*. 

Amplification  from  stimulated  emission,  offset  by  the  in- 
tracavity absorption  and  output  coupling  losses. 

• Formulation  in  terms  of  intracavity  photon  density: 

nph  ’ >/cl™ 


dn^lU/dt  ■ -ath]„phu)*^|-] 

• Cavity  buildup  time:  tc  • 


Fig.  2.4:  Basic  formulation  of  radiative  extraction  in  terms 
of  intracavity  photon  number  density,  with  buildup 
of  optical  fields  from  spontaneous  emission. 
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corresponding  to  the  laser  transition  from  an  excited  state  X to  a lower 
state  X,  as  well  as  by  absorption  processes  such  as 


(T _h  (A  ) * 

A + I(y)  abs  ,»  A 


can  be  written 

dl(t)/dt  = c[a(t)  - a^]  I(t  ),  (9) 

where  a(t)  is  the  instantaneous  net  gain  coefficient  in  the  medium  (i.e. 
active  gain  offset  by  absorption), 

a(t)  = [N*(t)  - N(t)]cr-  ^ nA(t)  °'abs(A)'  {10) 

A 

and  a is  the  cavity  threshold  loss  coefficient.  Actually,  (9)  applies 

only  to  the  situation  where  the  cavity  length  L and  the  gain  gength  L 

c 8 

of  the  active  medium  are  equal;  if  they  are  unequal,  it  must  be  modified 


to  read 


dl(t)/dt  = (L  /L  )c[a(t)-  ajKt) 
S ^ 


The  factor  L /L  has  the  significance  of  a contraction  in  the  radiative 
g c 

time  scale  to  account  for  the  fact  that  during  a round  trip  cavity  time, 
the  gain  medium  is  encountered  during  only  a fraction  (L^/L^)  of  that 
time.  Equivalently,  if  the  gain  medium  were  considered  to  be  distribu- 
ted uniformly  over  the  entire  cavity  length  L<c  (rather  than  over  the  ac- 
tual gain  length  L ),  then  the  actual  gain  would  have  to  be  reduced  by 

s 

the  factor  L^/L^  to  produce  an  equivalent  round  trip  amplification.  It 

should  also  be  noted  that  L /L  is  the  ratio  of  the  active  gain  volume 

g c 

to  the  optical  cavity  volume,  and  therefore  this  factor  logically  appears 
in  the  formulation  of  the  equation  for  conservation  of  power  density. 
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It  is  convenient  to  treat  stimulated  emission  and  absorption  processes 
in  a manner  parallel  to  that  for  kinetic  collision  reactions.  Therefore, 
the  equation  for  radiative  extraction  will  be  formulated  in  terms  of  the 
photon  number  density  in  the  cavity,  defined  by 

n = I/chv.  (12) 

ph 

This  will  also  facilitate  the  incorporation  of  spontaneous  emission  as  a 
source  term  for  the  buildup  of  the  optical  fields  from  noise.  The  equa- 
tion for  the  photon  number  density  includes  the  source  term  from  spon1- 
taneous  emission,  amplification  from  stimulated  emission,  and  losses 
from  output  coupling  and/or  intracavity  absorption: 

dnph/dt  = (Lg/Lc)  [ N*/t  (fi  / 4v)  + c[a(t)  - a^] nph]  . (13) 

a = [ y + (1 /2)  ln(  1 /R)  ] /L  is  the  threshold  loss  coefficient  expressed 
tn  g 

in  terms  of  the  intracavity  loss  per  pass  y and  the  output  coupling  re- 
flectivity R,  (n/4ir)  is  the  fractional  part  of  the  total  spontaneous  emis- 
sion which  builds  up  the  laser  mode  (taken  to  be  the  solid  angle  H sub- 
tended by  the  output  coupling  mirror  viewed  from  the  opposite  end  of  the 
cavity),  and  at  (t)  is  the  instantaneous  net  gain  coefficient  in  the  medium 
defined  by  Eq.  (10). 

Note  that  transient  effects  are  included  in  the  present  description  of  the 
radiation  field.  Although  a steady  state  oscillation  condition  (i.e.,  for 
which  the  gain  and  loss  coefficients  are  equal)  will  typically  occur  dur- 
ing the  time  scales  of  interest,  such  a condition  is  not  an  a priori  as- 
sumption of  the  present  analysis.  This  approach  allows  the  radiation  in- 
tensity to  be  described  for  situations  where  pumping  and  kinetic  times 
are  comparable  to  the  cavity  decay  time. 
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•c  * V w-  (14> 

For  such  situations,  gain  relaxation  effects  can  be  observed  in  which  the 
net  gain  initially  overshoots  the  cavity  threshold  and  subsequently  decays 
to  the  steady  state  gain  condition,  aft)  = a . Under  conditions  of  high 
pumping  the  overshoot  can  be  significant,  since  the  stimulated  emission 
from  the  medium  occurs  under  conditions  of  high  gain. 

When  formulated  in  terms  of  the  intracavity  photon  density,  the  radia- 
tion equations  can  be  included  in  a unified  way  with  the  other  kinetic  col- 
lision processes  of  the  laser  reaction  scheme.  For  example,  the  stimu- 
lated emission  process  (7)  can  be  described  by  the  equations 

dN*/dt  = - k.N*n  + k N n , = k(N-N*)nu  (15.  1) 

t pn  r pn  ph 

dN/dt  = - dN*/dt  (15.2) 

dn  , /dt  = (L  /L  )k(N*  - N)n  . , (15.3) 

ph  g c ph 

where  k^  = k^  = k = ccr,  and  cr  = stimulated  emission  (and  absorption) 
cross  section  (cm^).  Likewise,  the  contribution  of  reaction  (8)  to  the 
photon  field  is 

d>/dt=  ' (YLc,kabs<A>"A(t)V  <I6) 

where  k^^A)  = co‘aj3S(A)  is  the  rate  constant  for  absorption  (cm^/s). 

Although  minor  modifications  are  required  to  incorporate  the  factors 

(fi/47r)  and  (L  /L  ) and  the  term  - c a.,  n in  Eq.(13),  the  similarity 
g c tn  pn 

in  structure  of  Eq.  (3),  (15),  and  (16)  permits  the  radiation  field  (pho- 
tons) to  be  treated  on  an  equal  basis  with  the  other  molecular  species  in 
the  construction  of  the  master  equation  described  earlier. 
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2.3  Electron  Kinetics:  Boltzmann  Equation 


The  present  analysis  of  an  electrically  excited  laser  medium  incorpo- 
rates a completely  coupled  treatment  of  molecular  and  electron  kine- 
tics. The  electron  kinetics  analysis  is  based  upon  the  assumption  that 
the  plasma  has  attained  a quasisteady  state  on  the  time  scales  of  elec- 
trical pumping  and  inelastic  secondary  electron  collision  processes, 
so  that  the  electron  energy  distribution  function  f(u)  can  be  obtained 
from  solution  of  the  time- independent  Boltzmann  equation  as  a function 
of  gas  mixture  and  the  instantaneous  values  of  E/N  and  excited  state 
population  densities.  The  present  analysis  includes  electron- molecule 
inelastic  collisions,  elastic  electron- molecule  momentum  traisfer  col- 
lisions (with  recoil),  electron- electron  (i.  e. , Coulomb)  scattering,  and 
excitation  from  an  applied  (dc)  electric  field  as  well  as  source  terms 
for  electron  creation  (and  resulting  energy  deposition)  from  an  external 
e -beam.  Superelastic  electron  collisions  (i.  e.,  reverse  processes  in 
which  electrons  gain  energy  by  collisions  with  excited  molecular  spe- 
cies) are  important  in  situations  with  high  excited  state  population  den- 
sities, and  are  included  in  both  the  molecular  and  plasma  kinetic  equa- 
tions. Electron- electron  scattering  effects,  which  can  become  impor- 
tant under  conditions  of  high  fractional  ionization,  are  included.  From 
the  solution  of  the  Boltzmann  equation,  all  of  the  electron  transport  co- 
efficients, elastic  and  inelastic  secondary  electron  collision  rates,  el- 
ectrical power  partitioning,  etc.,  can  be  obtained. 

One  of  the  unique  features  of  the  present  analysis  is  a self-consistent 
treatment  of  secondary  ionization  and  attachment  phenomena  by  reten- 
tion of  a quasisteady  term  proportional  to  dne/dt  to  properly  account 
for  transient  plasma  effects  that  occur  on  the  time  scale  of  the  molecu- 
lar and  electron  kinetic  processes.  The  usual  formulation  of  the  time- 


ELECTRON  KINETICS  ANALYSIS 


Electron  Energy  Distribution  f0(u):  Numerical  solution  of  the  Boltzmann 
transport  equation. 

Quasisteady  state  approximation:  the  term  proportional  to  dne/dt  is  retained 
to  describe  situations  for  which  there  is  a net  creation  (or  loss)  of  electrons. 

Source  term  included  to  describe  external  ionization  (e.g.,  electron  beam  or 
photoionization): 

S(u)  ■ S0S(u)  + S(j(u). 

S0  term  provides  for  the  external  creation  of  electrons  at  zero  energy. 
Momentum  transfer  with  recoil;  elastic  heating  of  the  molecular  gas. 

Binary  inelastic  electron -molecule  collision  processes  included;  e.g., 

Vibrational  Excitation,  Electronic  Excitation,  Secondary 
Ionization,  Attachment,  Recombination,  Dissociation,  etc. 


• Electrons  created  by  secondary  ionization  assumed  to  be  at  zero  energy  and 
included  in  a source  term  ~V,-S(u). 

• Electron-electron  collisions  (Coulomb  scattering)  are  included  (important 
for  high  fractional  ionization,  £ 10'5). 

• Charged  particle  interactions  between  electrons  and  heavy  ions  neglected. 


• Superelastic  collisions,  important  for  excited  gas  mixtures,  are  included. 


• Electrical  parameters,  calculated  as  a function  of  E/N,  gas  mixture,  and  the 
excited  state  population  densities: 

Mobility,  Drift  Velocity,  Average  Energy,  Effective  Temperature, 
Characteristic  Energy,  Diffusion  Coefficient,  Plasma  Conducti- 
vity, Discharge  Current  Density,  Power  Densities. 

Forward  (and  reverse)  secondary  electron  excitation  rates. 

Fractional  partitioning  of  electrical  input  power  into  rate  of  energy 
storage  and  all  elastic  and  inelastic  scattering  mechanisms. 


Fig.  2.  5:  Summary  of  features  of  the  electron  kinetics  analysis. 
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independent  Boltzmann  equation  assumes  that  the  inelastic  processes 
include  only  binary  electron- molecule  collisions  for  which  there  is  no 
net  creation  or  loss  of  secondary  electrons.  When  creation  or  loss  of 
electrons  occurs  by  secondary  ionization,  attachment,  recombination, 
external  sources  of  ionization,  etc.,  it  is  necessary  to  retain  the  term 
proportional  to  dng/dt  in  the  time- dependent  Boltzmann  equation  in  or- 
der to  properly  describe  the  conservation  of  electron  density.  If  the 
electron  energy  distribution  f(u)  were  known,  the  rates  of  secondary 
ionization,  attachment,  recombination,  etc.  could  all  be  calculated  and 
dne/dt  could  be  determined.  However,  dne/dt  occurs  as  a parameter 
in  the  Boltzmann  equation,  and  in  order  to  calculate  f(u),  its  value  must 
be  known.  Therefore,  an  iterative  approach  was  developed  for  a self- 
consistent  determination  of  dne/dt.  The  numerical  techniques  for  solu- 
tion of  the  Boltzmann  equation  will  be  described  in  detail  later. 

The  basic  features  of  the  electron  kinetics  analysis  are  summarized  in 
Fig.  2.5.  The  scope  of  the  present  analysis  is  considered  to  be  quite 
comprehensive.  In  addition  to  the  fact  that  the  present  formulation  in- 
cludes more  refinements  than  many  of  the  Boltzmann  analyses  current- 
ly in  usage,  the  computer  subroutine  for  its  numerical  implementation 
is  quite  flexible  and  computationally  efficient. 

Theoretical  Formulation 

In  general,  the  electron  energy  distribution  function  for  the  gas  mix- 
tures and  relatively  low  E/N  values  typical  of  electric  discharge  lasers 
is  highly  non- Maxwellian.  The  increasing  availability  of  extensive  ex- 
perimental cross  section  data  for  elastic  and  inelastic  electron  collision 
processes  for  many  gases  of  interest  often  makes  it  possible  to  obtain 
a quantitative  analysis  of  electron  kinetics  by  direct  numerical  solution 
of  the  Boltzmann  equation.  If  there  is  no  net  creation  or  loss  of  elec- 
trons by  external  source,  secondary  ionization,  recombination,  or  at- 
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tachment,  the  electron  velocity  distribution  f(v)  for  a multicomponent 
gas  in  an  external  electric  field  E(r,  t)  satisfies  the  Boltzmann  trans- 
port equation, 

J^3/9t  +"v.v^  - ( e/m)  E (7,  t).  V-.  J f(v,  ?,  t)  = Sf/5t  (17) 

where  the  collision  term  on  the  right  hand  side  includes  the  effects  of 
all  elastic  and  inelastic  scattering  processes.  A detailed  discussion 

of  the  Boltzmann  transport  equation  can  be  found  in  works  by  Holstein^ 

8 9 10  11 
Allis  , Schkarofsky  ^t  al  , Frost  and  Phelps  , Carleton  and  Megill  , 

12  13 

and  Englehardt  et  al  ’ . Various  approximations  and  techniques  for 

its  numerical  solution  and  application  to  the  analysis  of  electric  dis- 

14  15 

charge  gas  lasers  have  been  discussed  by  Nighan  £t  al  , Rockwood  , 

16  17  18  19 

Hancock  et  al  , Elliott  £t  al  , Morgan  and  Fisher  , Lowke  £t  al  , 

Lacina^,  and  others.  The  present  approach  is  based  upon  extensions 

and  generalizations  of  earlier  work^  related  to  analysis  of  CO  kinetics. 

At  the  gas  densities  of  interest  for  laser  operation,  elastic  and  inelas- 
tic electron- molecule  collision  frequencies  are  so  high,  relative  to 
those  for  molecular  kinetic  processes,  that  the  electron  distribution 
function  can  be  assumed  to  have  adjusted  "instantaneously"  on  the  time 
scales  over  which  the  excited  state  population  densities  are  evolving. 

It  follows,  therefore,  that  the  first  term  of  Eq.  (17)  containing  the  par- 
tial time  derivative  can  be  discarded  if  the  electric  field  E(r,  t)  is  as- 
sumed to  be  dc  or  slowly  varying.  For  a spatially  uniform  field,  it  is 
similarly  possible  to  neglect  the  second  term,  since  the  mean  free 
path  (which  is  related  to  the  elastic  collision  cross  section  and  mole- 
cular density)  is  typically  small  compared  to  the  discharge  dimensions. 
Thus,  the  Boltzmann  equation  (17)  becomes 


- (e/m)  E*  V-f(v)  = S f / 6 1 J 
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The  usual  approach  for  solving  this  equation  begins  with  an  expansion  of 
f(v)  into  Legendre  polynomials, 


f(v)  = fQ(  v)  + ( v/v)  . fx  ( v)  + [3  vv/v2  - TJ:  f^(  v)  + . . . (19) 

( 1 and  f^  represent  dyadic  quantities.)  Usually,  f can  be  well  approxi- 
mated using  only  the  first  two  terms, 

f (v)  = fQ(v)  + (v/v).  f^v).  (20) 

Physically,  the  justification  of  a two- term  approximation  is  related  to 
assumptions  about  the  degree  of  anisotropy  to  be  expected  for  the  velo- 
city distribution.  If  the  electric  field  is  small  enough  that  the  directed 
speed  of  the  electrons  (as  measured  by  their  drift  velocity)  is  much  less 
than  their  random  thermal  velocities,  the  small  first  order  perturbation 
v . fj  from  an  isotropic  distribution  fQ  ( v)  should  be  a good  approxima- 
tion. For  high  values  of  E/N  characteristic  of  self- sustained  discharges 

21  ~ 

and  for  certain  gases,  it  has  been  speculated  that  the  f 2 term  may  be 
comparable  to  the  f j term,  although  no  analysis  has  yet  been  undertaken 
to  include  these  higher  order  terms  quantitatively.  Retention  of  more 
than  the  first  two  terms  would  require  knowledge  of  angular  (i.  e.  differ- 
ential scattering)  cross  section  data,  which  is  not  generally  available. 
Typical  experimental  data  gives  only  integrated  cross  sections  (over  all 
angles)  as  a function  of  the  electron  energy.  Some  of  the  available  data 
in  the  literature  has  been  determined  indirectly  by  fitting  measured 
transport  coefficients  to  synthetic  cross  sections,  using  a numerical 
Boltzmann  analysis  based  upon  the  two- term  expansion  (20)  of  f(v). 

Thus,  consistent  use  of  that  data  would  require  that  the  electron  kine- 
tics calculations  be  based  upon  the  two- term  approximation.  For  values 
of  E/N  typical  of  e -beam  sustained  discharges,  this  approach  is  justi- 
fied, and  it  is  numerically  much  simpler  than  would  be  the  case  if  angu- 
lar calculations  were  required. 
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In  the  expansion  (19),  fQ(v),  f^v),  ^(v),  . . . are  functions  only  of  the 
magnitude  v =|  v),  since  all  of  the  angular  dependence  has  been  put  ex- 
plicitly  into  the  spherical  harmonic  coefficients  v,  (3w  - 1),  ...  etc. 

It  can  easily  be  verified  that 


V f(v)  = — —■  (d/dv)(v2f  ) + (df  /dv)v  + (v/3)  (d/dv)(  v f ) • 
v 3 v4  l o 1 


(3 w - 1)  + ... 


(21) 


where  V^f(v)  has  been  separated  into  spherical  harmonic  terms  which 
have  l-  0,  1,  2 symmetry,  respectively.  If  (21)  is  substituted  into  (18) 
and  separated  into  scalar  and  vector  parts  (which  corresponds  to  multi- 
plying by  1 or  v and  integrating  over  the  spherical  solid  angle  fi^),  the 
following  equations  result: 


- (eE/3  mv^)  • (d/dv)  ( v2  f ) = 8f  /8t 

1 o c 


(eE/m)(df  /dv)  = S f /St 

o 1 c 


(22.1) 


(22.  2) 


The  most  important  physical  mechanism  for  reducing  the  asymmetry  in 

the  distribution  function  is  electron- molecule  momentum  transfer  colli- 

— 9 

sions.  The  collisional  rate  of  change  of  f^  can  be  approximated  by 


= - 2 


(23) 


X X 

where  V (v)  = NvvQ  (v)  is  the  momentum  transfer  collision  frequen- 
m A m 

cy  for  molecules  X in  terms  of  the  population  density  N and  momentum 

X X 

transfer  cross  section  Qm  (v).  In  terms  of  the  electron  energy 


u = mv  / 2e 


(24) 


(with  u in  units  of  eV,  e = 1. 602  10"  12  erg/eV  = 1. 602  10*  ^ J/eV),  the 
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equation  for  f (u)  becomes,  after  manipulation  of  Eq.  (22)  and  (23) 


df  (u)/du 


Expressions  for  the  collision  term  8f  /St  on  the  right  hand  side  of 

°7  c 8 9 

Eq.  (25)  have  been  derived  by  Holstein  , Allis  , Schkarofsky  et  al  , 

Frost  and  Phelps^,  and  numerous  others.  It  shall  not  be  the  purpose 
of  the  present  discussion  to  repeat  these  extensive  and  rigorous  deri- 
vations here. 


However,  it  may  be  useful  at  least  to  motivate  the  form  of  the  results 
for  inelastic  binary  collisions  characterized  by  an  isotropic  scattering 
cross  section.  For  an  inelastic  scattering  process,  the  collision  term 
is  defined  simply  as  the  net  flux  of  electrons  scattered  into,  and  out  of, 
a region  £v",  v + A"v}  of  velocity  space.  The  electrons  scattered  into 
[v,  "v  + Av]  originate  from  all  volume  elements  [ v’1,  v*  + Av1  ] for  which 
energy  is  conserved: 


where  u is  the  inelastic  energy  loss  (eV)  for  the  scattering  process, 
o 

and  Eq.  (24)  was  used.  Under  the  assumption  of  isotropic  scattering, 

Eq.  (26)  uniquely  defines  Av'  for  an  arbitrary  choice  of  v and  Av,  for 
2 ^ 2 

we  can  set  Av  ~ v Av,  Av1  ~ v'  Av1,  and  invoke  Eq.  (26)  to  obtain 
Av'/Av  = v/v1.  Therefore,  it  follows  that 


mv 

8£o 

2e 

8 t 

Av'  = (v'/v)a7. 


(27) 


The  rate  (cm  *)  at  which  electrons  are  scattered  out  of  [v'.v'+Av1] 
is  given  by  f(v')Av'N  v'Q(v'),  where  Q(v')  is  the  scattering  cross  sec- 
tion and  N is  the  neutral  molecule  particle  density.  Recall  that  f(v')Av' 
represents  the  spatial  density  of  electrons  (cm  in  the  velocity  region 
P.  7' + Av'],  and  that  Nv'Q  (v1)  has  the  significance  of  collision  fre- 
quency (s  In  order  to  obtain  an  expression  for  the  total  flux  of  elec- 

trons scattered  into  ["v,  7+A\T]  from  [v1,  v'+  Av1],  it  is  necessary  to 
sum  the  contributions  from  all  velocities  v'  which  satisfy  the  conserva- 
tion of  energy  condition.  The  required  "sum"  over  the  continuum  of 
velocities  v',  restricted  by  Eq.  (26),  can  be  expressed  as  an  integral 
over  dv',  weighted  by  a density  function 


F(u')  = 6 (u'  - u - u ) / f dv"  8 (u"  - u - u ). 

o J o 


It  is  easy  to  show  that  F(u')dv'  = du'  8(u'  - u - u ),  and  therefore  the 

o 

restricted  "sum"  over  velocities  v1  is  equivalent  to  integration  over 
all  energies  u'  weighted  by  the  distribution  6(u'  - u - uq).  Thus,  the 
collision  term  can  be  written 


Av  8f(v)/St  | = - f (v) Av  v N Q(v) 


j du'8(u'  - u-  uq)  f (v')Av' v1  NQ(v'), 


where  the  first  term  represents  scattering  out  of  (and  the  second  term 
scattering  into)  the  region  [v,  v + ATv  ] of  velocity  space.  From  Eq.  (27) 
it  follows  that  the  collision  term  can  be  expressed  as 
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8f(v)/8t  = (N/v)  [f  (v1)  v'2  Q(v')  - f (v)  v2  Q (v)l  (28) 

c v'  = v + 2eu  / m 


To  lowest  order  in  the  expansion  (20),  f(v)  = fQ(v)  depends  only  upon 
the  magnitude  v,  as  is  the  case  for  Q (v),  and  therefore  all  quantities 
can  be  expressed  as  functions  of  the  energy  u.  Thus,  the  lowest  order 
inelastic  collision  term  becomes 


Sf  (u)/8t  = — 

o mv 

c 


N f(u+u  )f  (u  + u )Q(u  + u ) - uf  (u)Q(u)l,  (29) 

L o o o o o J 


where  Eq.  (24)  and  (28)  have  been  used.  For  all  collision  processes 

(both  elastic  and  inelastic),  there  will  be  a contribution  to  8f  / S t I 

o | c 

originating  from  momentum  transfer.  (Whenever  an  electron  of  energy 
u undergoes  a collision  process  with  a neutral  molecule  of  mass  M, 
there  is  a transfer  of  energy  of  order  (m/M)  u resulting  from  mo- 
mentum conservation.)  The  momentum  transfer  contribution  is^ 


8f  (u)/8t  = — 

o mv 

c 


<*» 


Combining  Eq.  (25),  (29),  and  (30),  we  obtain  the  following  form  for  the 
Boltzmann  equation: 

(E2/3)  d/du£udfQ/du/<  NQm>]  +d/du  £u2<(2m/M)NQm>  (fQ  + ^ -j~)  j 

+ ^[(u  + u«)£0(u  + u«)NaQa^u  + ua)  " ufQ(u)  N«Qo((u) 
a 

+ (u  - ua)fo(u  -ua)N*  Q_a(u-ua)  - ufo(u)N*  Q.a(u)]  = o (31) 
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where  Q±  a(u)  are  the  (forward  and  reverse)  cross  sections  for  the  in- 
elastic  collision  processes,  labeled  by  a , Na  and  Na  are  the  population 
densities  of  the  lower  and  upper  (excited)  states,  and  ua  is  the  inelastic 
energy  loss  associated  with  the  collision.  The  two  weighted  averages 
which  occur  in  Eq.  (31)  are  defined  by 

< N Q ( u) > = V N Q*(u),  (32.1) 

m t—i  A.  m 

X - 

<(2m/M)NQm(u)>  = £ (2m/Mx)  Nx  Q^(u)  (32.2) 

X 

The  principle  of  detailed  balance  provides  relations  between  the  cross 
sections  Q±<j(u)  for  forward  and  reverse  ( super  elastic)  processes: 

(u  + ua)Qa(u  + ua)  = uQ_a(u)  (33.1) 

(u  - ua)Q_a(u  - ua)  = uQa(u)  (33.2) 

The  first  two  terms  in  the  sum  over  a in  Eq.  (31)  correspond  to  colli- 
sions in  which  the  electrons  lose  energy,  while  the  last  two  terms  (which 
represent  the  superelastic  processes)  correspond  to  collisions  in  which 
excited  molecular  species  transfer  their  energy  to  the  electrons.  Only 
binary  electron- molecule  collisions,  including  vibrational  excitation, 
electronic  state  excitation,  attachment,  recombination,  ionization,  etc. 

will  be  considered.  It  should  be  emphasized,  however,  that  the  formu- 
lation of  the  Boltzmann  equation  as  given  by  Eq.  (31)  applies  only  to  a 
quasisteady  state  situation  for  which  there  is  no  net  creation  or  loss  of 
electrons.  If  processes  such  as  secondary  ionization,  electron  attach- 
ment or  recombination  are  included,  Eq.  (31)  is  valid  only  under  condi- 
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tions  for  which  the  collision  frequencies  for  such  processes  are  much 
lower  than  those  for  the  inelastic  excitation  processes  which  do  not  in- 
volve the  creation  or  loss  of  electrons.  It  should  be  noted  that  the  col- 
lision term  for  such  inelastic  processes  does  not  contain  all  four  of  the 
terms  indicated  in  the  general  expression  given  in  Eq.  (31).  For  ex- 
ample, only  the  second  term  occurs  for  recombination  or  attachment 
(in  which  an  electron  is  lost),  and  for  ionization  or  dissociation,  only 
the  first  two  terms  occur  (since  there  is  no  reverse  binary  collision 
process).  Because  of  the  low  density  of  states  available  for  momentum 
conservation,  two-body  recombination  is  generally  a very  slow  process. 
Electron  recombination  proceeds  as  a three-body  process,  which  can- 
not be  treated  according  to  the  present  formulation.  However,  it  is  of- 
ten found  experimentally  that  the  rates  for  three-body  recombination 
vary  with  pressure  as  if  the  process  were  two-body,  and  thus,  effective 
cross  sections  can  be  used  to  describe  it  as  a binary  collision  appropri- 
ate for  the  formulation  of  Eq.  (31).  Thus,  for  (schematic)  processes 
such  as 

Attachment:  A + e — *-  B + C (34.  1) 

Two- Body  Recombination:  A+  + e — *•  B + C (34.2) 

the  inelastic  collision  term  in  Eq.  (31)  consists  only  of 

- ufo(u)NaQa(u).  (35) 

For  processes  such  as 

Secondary  Ionization:  B + e — ► B*  + e + e (36.  1) 

Dissociation:  AB  + e — ► A + B + e (36.  2) 

the  collision  term  contains  only 
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(u  + via)fo(u  + ua)NaQa(u  + ua)  - ufQ(u)Na  Qa(u).  (37) 

Furthermore,  in  the  case  of  secondary  ionization,  there  is  the  crea- 
tion of  an  additional  electron.  Because  experimental  knowledge  of  io- 
nization cross  sections  is  not  usually  available  in  the  detail  necessary 
to  describe  the  final  energy  distribution  of  the  created  electrons,  it  is 
often  assumed  that  the  additional  electrons  are  created  at  zero  energy, 
and  that  they  relax  into  the  distribution  on  the  same  time  scale  (given 
by  the  elastic  collision  frequency)  over  which  the  rapid  attainment  of 
a quasisteady  state  occurs.  Therefore,  the  collision  term  for  secon- 
dary ionization  requires  an  additional  modification,  the  inclusion  of  a 
term 

co 

8(u)Na  j*  duufQ(u)  Qa(u).  (38) 

o 

Perhaps  a more  realistic  hypothesis  would  be  to  assume  that  the  final 
electrons  have  equal  energies  or,  better  still,  that  their  energies  are 
partitioned  according  to  a purely  kinematic  distribution  of  a three- 
body  breakup.  It  is  believed  that  such  assumptions,  which  needlessly 
complicate  the  analysis,  would  provide  no  additional  physical  insight. 
Therefore,  the  electrons  created  by  secondary  ionization  are  assumed 
to  have  zero  energy. 

The  effects  of  recoil  in  momentum  transfer  collisions  are  included  in 
the  second  term  of  Eq.  (31),  where  the  molecular  energy  distribution 
has  been  assumed  to  be  Maxwellian  at  temperature  T.  Note  that,  if 
there  is  no  electric  field  and  if  kT  is  much  less  than  the  energy  thresh- 
olds for  the  inelastic  collision  processes,  then  an  approximate  solu- 
tion for  the  electron  energy  distribution  is  obtained  by  equating  the 
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second  <-erm  of  Eq.  (31)  to  zero,  giving 


K 


. 


fQ(u)  = exp  [-eu/kT]. 


This  is  to  be  expected,  since  electrons  and  molecules  would  be  in  ther- 
mal equilibrium  with  each  other  in  that  case.  The  momentum  transfer 
term  in  Eq.  (31)  is  important  for  low  values  of  E/N  at  low  tempera- 
tures, and  for  gas  mixtures  containing  light  molecules.  The  physical 
significance  of  that  term  is  elastic  heating  of  the  molecular  gas,  to  be 
discussed  in  more  detail  in  Sec.  2.  5.  The  higher  the  value  of  the  low- 
est inelastic  excitation  threshold,  the  more  important  will  be  the  frac- 
tional power  partition  into  elastic  heating  for  a given  value  of  E/N. 

For  molecular  gases,  characterized  by  vibrational  level  excitation 
processes  with  very  low  energy  thresholds  (typically  ^ 0.  1 - 0.  2 eV), 

elastic  heating  is  negligible  for  values  of  E/N  characteristic  of  typical 

- 17  2 

laser  excitation  ( J2  3 x 10  Vcm  ).  However,  for  rare  gas  laser  mix- 
tures for  which  the  lowest  inelastic  energy  loss  usually  corresponds  to 
excitation  of  an  electronic  level  (~  few  eV),  the  contribution  of  elastic 
heating  can  be  significant  even  for  relatively  high  values  of  E/N  (e.  g. 

~ 10'^Vcm^).  Effects  of  rotational  excitation,  which  are  not  in- 
cluded in  the  present  version  of  the  electron  kinetics  analysis,  can  be 
9 10 

incorporated  * into  the  Boltzmann  equ  vion  in  the  same  way  as  mo- 
mentum transfer  by  making  a suitable  continuum  approximation. 

Most  authors  have  not  only  restricted  their  attention  to  solution  of  the 
Boltzmann  equation  in  the  form  of  Eq.  (31),  but  have  often  ignored  the 
superelastic  collision  terms.  For  weakly  excited  gases,  such  an  ap- 
proach is  justified,  providing  that  attachment  and  secondary  ionization 
could  be  neglected.  However,  in  order  to  properly  describe  effects  of 
interest  in  electrically  excited  lasers,  superelastic  collisions  must  be 


I 
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included,  and  in  addition,  several  important  extensions  to  the  formula- 
tion (31)  of  the  Boltzmann  equation  are  necessary. 

First  of  all,  the  contribution  (38),  representing  the  creation  of  electrons 
(at  zero  energy)  from  secondary  ionization,  must  be  included  in  the  col- 
lision term  of  Eq.  (31). 

Secondly,  it  is  useful  to  incorporate  a source  term  Sgxt(u)  to  describe 
(as  a function  of  energy)  the  creation  of  electrons  by  external  sources  of 
ionization  (e.  g.  , e -beam,  photoionization,  etc.  ) and  the  resulting  ener- 
gy deposition.  Since  there  are  often  several  electron  creation  processes 
for  which  no  a priori  knowledge  of  this  energy  dependence  is  available 
(e.  g.  , Penning  ionization,  photodetachment,  photoionization  by  intraca- 
vity radiation),  it  is  convenient  to  assume  that  the  resulting  electron  en- 
ergy is  zero.  Thus,  the  source  term  is  assumed  to  be  of  the  form 


S (u)  = S 8(u)  + S (u),  (39) 

ext  o 

-3-1 

where  S represents  a rate  (cm  s ) of  creation  of  electrons  at  zero 

° -3-1 

energy,  and  S(u)  represents  a rate  of  creation  per  unit  energy  (cm  s 

eV  *)  of  electrons  of  energy  u (i.  e.  , S(u)  du  is  the  rate  of  creation  of 
electrons  in  the  energy  interval  [u,  u + du]).  The  first  term,  which  can 
be  used  to  include  all  electron  creation  processes  arising  from  excited- 
state  interactions  in  the  medium,  contributes  nothing  to  energy  deposi- 
tion into  the  electron  gas  (since  the  electrons  have  zero  energy).  The 
second  term  S (u)  is  used  to  describe  electron  creation  from  an  exter- 
nal source,  typically  e’-beam  ionization.  A more  detailed  discussion  of 
energy  (power)  balance  is  given  in  Sec.  2.  5. 
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Third,  in  order  to  properly  conserve  electron  number  density,  a term 
proportional  to  dne/dt  must  be  retained  from  the  time- dependent  Boltz 
mann  equation.  The  term  dne(u,  t)/dt  (where  ne(u,  t)  du  =ne(t)duu 
fQ(u)  is  the  electron  density  in  the  energy  interval  [u,  u + du])  can  be 
written 


dn^fu,  t)/dt 


1 / 2 

u fQ(u,  t)  dne/dt 


+ n^tju1^2  dfQ(u,  t)/dt. 


(40) 


It  is  reasonable  to  retain  the  earlier  assumption  that  the  shape  fG(u,  t) 
is  not  an  explicit  function  of  the  time,  since  a quasisteady  state  will  be 
established  on  a time  scale  much  shorter  than  that  characteristic  of  mo- 
lecular and  electron  kinetic  processes.  Thus,  fQ(u)  will  be  only  an  im- 
plicit function  of  time  through  its  dependence  upon  instantaneous  values 
of  various  parameters,  such  as  E/N,  and  the  excited  state  population 
densities  Na(t).  It  follows  that  the  second  term  in  Eq.  (40)  can  be  de- 
leted, and  the  explicit  time  dependence  in  fQ(u, t)  suppressed  in  the  first 
term.  The  first  term  must  be  retained,  even  in  the  quasisteady  state 
approximation,  since  it  represents  the  net  rate  of  creation  (or  loss)  of 
electrons  on  the  time  scale  of  attachment,  recombination,  ionization, 
etc.  It  will  be  shown  that  retention  of  the  dne/dt  term  is  necessary  to 
insure  that  the  first  integral  of  the  Boltzmann  equation  correctly  results 
in  the  conservation  equation  for  the  electron  number  density  nft(t). 

Since  dne/dt  occurs  as  a parameter  whose  value  cannot  be  determined 
until  the  equation  is  solved,  a self-consistent  iterative  method  (to  be 
discussed  in  more  detail  later)  was  developed. 


Finally,  the  effects  of  Coulomb  scattering  can  become  important 
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under  conditions  of  high  fractional  ionization,  so  electron- electron  col- 
lisions must  be  included.  (Charged  particle  interactions  between  elec- 
trons and  heavy  ions  will  be  neglected  in  the  present  analysis.  ) For  a 
highly  ionized  gas,  electron- electron  collisions  dominate  electron- mole- 
cule collisions,  and  this  tends  to  force  the  electron  energy  distribution 
to  become  more  nearly  Maxwellian.  Rockwood^  has  shown  how  the 
Boltzmann  equation  can  be  expressed  in  terms  of  a "flux  divergent"  de- 
scription of  the  electron  distribution  in  energy  space,  and  how  electron- 
electron  collisions  can  be  incorporated  into  such  a formulation. 


Accounting  for  these  refinements,  the  Boltzmann  equation  becomes 


u ' f (u)  dn  /dt  + dj  (u)  / du  + dj  ,(u)/du  + dj  (u)/du  = 

O C t 61  6“ 

00 

SQ  8 (u)  + S (u)  + ne |s(u)£Ni  Jduu  fQ(u)  Q.  (u)  + 

^ i o 

^ j^(u  + ua)fo(u+ua)NttQa(u  + ua)  - ufQ(u)Na  Qa(u) 
a 

+ (u-  ua)fQ(u-  ua)N*  Q a (U-  Ua  ) - ufQ(u)N*  Q ^ (u)]  J . 


The  sum  over  i represents  a sum  over  all  secondary  ionization  proces- 
ses. The  sum  over  or  corresponds  to  the  complete  set  of  all  inelastic 
electron- molecule  collisions,  including  ionization,  with  the  tacit  under- 
standing (according  to  the  foregoing  discussion)  that  certain  terms  in 
the  general  collision  expression  are  to  be  omitted  for  processes  for 
which  there  is  a net  creation  or  loss  of  electrons,  or  for  which  there  is 
no  reverse  binary  collision.  (I.  e.  , the  collision  term  is  reduced  to  the 
expressions  (35)  or  (37)  in  such  cases. ) The  "current  densities"  Jf(u), 
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Jel(u),  and  Jee(u),  which  correspond  to  a flux  of  electrons  in  energy- 
space  driven  by  the  applied  electric  field,  elastic  scattering  collisions 
and  electron- electron  collisions,  are  defined  by 


Jf(u)  = - (EZ/3)ne>JS  [ udfo/du/<NQm(u)  >]  (42.1) 

“!<l2»'M|NQm(u|  > [fQ  4 ] <42.  2, 

J (u)  = - (2W3)n2\/|?(q4/e2)lnA[P(u)£°  + Q(u)  f ] (42.3) 

ee  c V m du  ° 

where 

u ao 

P(u)  = 2 f dww^2fQ(w)  + 2 u2^2  J*  dw  fQ(w), 
l o u 

U 

Q (u)  = 3 J*  dww'  2fQ(w). 
o 

- 10 

In  the  above  equations,  q = 300  e = 4.8  x 10  esu  is  the  electronic 

- 12 

charge  (cgs  units),  e * 1.  602  x 10  erg/e V,  and  A is  the  ratio  of  the 
Debye  length  and  the  classical  distance  of  closest  approach  rQ  (for 
an  electron  of  average  energy  u = u ),  defined  by 


(43. 1) 

(43.2) 


A = D/ro  (44) 

lD  = (kTe/4irn*)1/2  (45) 

r = q2 /( eu ) = 2q2/3kT  . (46) 


Comparison  of  Eq.  (31)  with  Eq.  (41)  shows  that  they  are  equivalent, 
except  for  the  inclusion  of  the  additional  terms  dne/dt  for  electron  con- 
servation, S(u)  for  external  ionization  sources,  dJee(u)/du  for  electron- 
electron  collisions,  and  the  term~8(u)  representing  the  creation  of  el- 
ectrons (at  zero  energy)  by  secondary  ionization  and  the  external  source 
SQ.  The  formulation  (41)  of  the  Boltzmann  equation  to  be  used  in  the  pre- 
sent analysis  is  written  in  such  a way  that  its  significance  as  a continuity 
equation  (in  energy  space)  is  manifestly  apparent.  In  Fig.  2.6,  an  in- 
terval [u,  u + du]  of  energy  space  is  depicted,  with  an  electron  "current 
density"  J(u).  If  there  is  an  external  source  S ext(u)  and  a (non-local) 
source  SCQ^(u)  from  inelastic  collisions  driving  electrons  into  (and  out 
of)  the  interval  [u,  u + du],  the  continuity  equation  can  be  written 

du[dne(u,  t)/dt]  + [ J(u  + du)  - J(u)]  = du  [Sext  (u)  + Scou(u)] 
or, 

dng(u,  t)/dt  + dJ(u)/du  = StQt(u)  (47) 


J (u)  J (u  + du) 


u u + du 


Fig.  2.  6:  Continuity  equation  interval. 

Thus,  the  first  term  on  the  left  hand  side  of  Eq.  (41)  represents  the 
rate  of  change  of  electron  density  in  the  interval  [u,  u + du],  while  the 
dJ/du  terms  represent  electron  flux  in  energy  space  driven  by  the  ap- 
plied electric  field,  elastic  collisions,  and  electron- electron  scatter- 
ing. On  the  right  hand  side  of  (41),  there  are  source  terms  from  ex- 
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ternal  ionization,  secondary  ionization,  and  nonlocal  inelastic  scatter- 
ing of  electrons  into  and  out  of  the  interval  [u,  u + du]. 


The  normalization  condition  for  f (u)  is  taken  to  be 

o 


oo 


f 


du  u1/2  f (u)  = 1, 

o 


(48) 


1/2 


where  the  factor  u originates  from  the  density  of  states  associated 
with  the  transformation  from  the  three-dimensional  velocity  spaced  to 
the  energy  space  u:  d v = 47rv^dv  = 27r(2e/m)^  ^ u^^  du.  (Thus,  the 
units  of  fQ(u)  will  be  eV"^^. ) From  the  solution  fQ(u)  of  the  Boltzmann 
equation  (41),  all  of  the  forward  and  reverse  electron  excitation  rates 
<vQ±a(u)>,  electrical  power  partitioning,  and  miscellaneous  plasma 
parameters  (drift  velocity  v^,  mobility  fj.,  average  and  characteristic 
energies  u and  €^»  effective  temperature  Tg,  etc.)  can  be  obtained: 


oo 


Collision  Rates  (cm  /s): 

<CvQ±a('i)  > = (2e/m)1^2  J du  u fQ(u)Q±a(u)  (49) 

o 

Mobility  (cm2/Vs): 

oo 

fjL  = - (l/3)(2e/m)1/2  f du[u/<NQ  (u)>]  df  /du  (50) 

J m ' o 


Diffusion  Coefficient  (cm  /s): 


oo 


D = (1/3)  (2e/m)^2  fdufu/^NQ  (u)S]  f (u)  (51) 

J m o 

o 

Average  Energy  (eV): 


oo 


u = J"  du  u^2fQ(u) 


(52) 
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Drift  Velocity  (cm/s) 


For  the  case  of  an  ac  electric  and/or  applied  magnetic  fields,  there 

are  other  transport  coefficients  of  interest,  and  expressions  for  these 

quantities  as  well  as  discussion  of  the  extensions  of  the  Boltzmann 

equation  required  to  describe  these  phenomena  can  be  found  in  the  ref- 
8-10 

erences 


The  first  integral  (from  0 to  oo)  of  Eq.  (41)  gives  an  expression  for 
the  conservation  of  electron  density: 


are  the  total  secondary  ionization  and  electron  attachment  frequencies 
Note  that  for  all  of  the  inelastic  scattering  processes  which  conserve 
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electrons,  the  integral  over  the  collision  term  vanishes: 
co 

Jdu  [(u+ua)fQ(u+aa)Qa(u  + ua)  - ufo(u)Qa(u)]  = Q,  (59) 
o 

which  can  be  shown  by  shifting  the  variable  of  integration  and  noting 
that  Qa(u)  = 0 for  u<ua.  For  attachment  or  recombination,  however, 
only  the  term  (35)  occurs  in  the  collision  expression,  leaving  the  v 

2L 

term  in  Eq.  (56).  For  secondary  ionization,  both  terms  shown  in  Eq. 
(59)  appear  and  the  integral  vanishes:  however,  there  is  still  a term 
in  Eq.  (56)  containing  originating  from  the  5(u)  term  in  Eq.  (41). 
The  integral  of  S(u)  on  the  right  hand  side  represents  the  total  rate  of 
change  of  electron  density  from  the  external  source. 

It  is  apparent  that  the  term  proportional  to  dne/dt  must  be  retained  in 
order  to  correctly  obtain  conservation  of  electron  density  in  situations 
where  creation  or  loss  of  electrons  is  important.  Obviously,  this  will 
always  be  the  case  for  transient  analysis  of  pulsed  electrical  lasers  if 
it  is  important  to  properly  describe  the  build-up  of  the  plasma.  Note 
that  Eq.  (31),  which  does  not  include  the  dne/dt  term,  may  give  rise 
to  solutions  for  fQ(u)  which  contain  logarithmic  singularities  at  u = 0. 
To  demonstrate  this,  consider  the  case  of  attachment  only,  and  inte- 
grate the  Boltzmann  equation  from  u to  oo: 

J (oo)  - J(u)  = - n v 

6 d 

As  u-»0,  (since  J (oo)  = 0),  this  becomes  an  equation  of  the  form 

[Au^  + Buldf  /du  = n v 
L J o e si 

which  is  characterized  by  a singularity  at  u = 0. 
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If  Eq.  (41)  is  multiplied  by  eu.  and  integrated  over  all  energy,  an  ex- 
pression for  electrical  power  balance  is  obtained: 


VJU 

n^/iE2  + j*du(eu)S(u)  = e u dne/dt  + 


ne  e(2e/m)1/2  J J du  u2<(2m/M)  N [fQ  + (kT/e)dfQ/du  j 
o 

00  ) 

+ 2ua  J*du  u fQ(u)  [NaQ«(u)  - Na  Q_a(u>]  j (60) 

The  form  of  the  inelastic  collision  term  on  the  RHS  of  Eq.  (60)  results 
from  a change  of  variables  in  the  first  term  of  an  expression,  such  as 

(Co  °°  ) 

No|J*du  u(u+  ua ) fQ(u+  ua ) Qa(u  + ua)  - J du  u2fQ(u)  Qa(u) 

o o f 

oo 

= Na  Jduf^u)  Qa(u)  C u(u  • ua)  “ u2  ] 

o 

oo 

= -NauaJ*  du  u fo(u)Qa(u)  (61)  . 

o 

Note  that,  for  processes  such  as  ("two- body")  electron  recombination 
and  attachment  (in  which  an  electron  is  lost),  only  the  second  term  of 
the  form  (35)  (i.e.,  - ufQ(u)  Na  Qa(u))  occurs  in  Eq.  (41).  In  that 
case,  the  first  term  on  the  LHS  of  Eq.  (61)  is  absent  and  the  corres- 
ponding term  in  the  power  balance  equation  (60)  must  be  modified: 

oo 

Na  Jdu  u2fQ(u)  Qa(u)  (62) 

o 

The  physical  interpretation  of  Eq.  (60)  is  that  the  total  electrical  input 
power  density  (discharge  and  external  deposition)  on  the  left  hand  side 
is  partitioned  into  eudne/dt  (representing  the  rate  of  change  of  electron 
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kinetic  energy  stored  in  the  electron  gas),  elastic  heating  from  moraen- 
turn  transfer  collisions,  and  dissipation  by  the  inelastic  excitation  pro- 
cesses on  the  right  hand  side.  Note  that  the  contribution  from  electron- 
electron  scattering  vanishes- -this  mechanism  merely  redistributes  the 
electron  energy  and  tends  to  drive  the  distribution  toward  a Maxwellian 
shape.  To  show  that  the  electron- electron  contribution  vanishes,  inte- 
grate by  parts, 


oo 


oo 


oo 


'/duJee(u), 


(63) 


/duu  dJ  (u)/du  = J (u) 
ee  ee 

O O 

and  the  first  term  on  the  right  vanishes  because  P(oo)  = 2 u,  Q(oo)  = 3, 

P(0)  = Q(0)  = 0,  and  f (oo)  = df  (co)/du  = 0.  Thus, 

o o 

co  oo 

- J*du  u dJee(u)/du  = J*du  J (u) 


ee 


oo 


(2ir/3)  n^(q4/e2)  InA  J*  du[PdfQ/du  + Q fQ(u)]  (2e/m) 1 


Suppressing  the  constant  factor,  and  integrating  the  first  term  by  parts 
gives 


oo 


j oufPdfo/du  + Qf  ] = P(u)f  (u)  + fduf  [Q  - dP/du] 

o oo  ° ° o 

= j* du  fQ(u)  f Q(u)  - dP/du]  (64) 

o 

From  the  definitions  (43), 


u 


00 


Q(u)  - dP(u)/du  = 3;  J*dww^2fo(w)  - u1^2  j*dwfQ(w)  J (65) 


I 


and  the  expression  in  (64)  becomes 
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00 

J du  fQ(u)  [Q(u)  - dP(u)/du]  = 


oo  u oo  oo 

3 j J*dufQ(u)  j*  dww  ^fQ(w)  - J*  duu^^fju)  J*  dwf^(w)  j 


Interchanging  the  variables  (u,  w)  in  the  second  term  above  gives 


oo 


j*  du  fQ(u)  [ Q (u)  - d P (u)  /du  ] = 


w 00  u oo  oo 

o o O W 

Since  both  of  the  integrals  in  Eq.  (66)  are  carried  out  over  the  same 
(infinite)  45°  sector  in  the  upper  half  of  the  (u,  w)  plane  defined  by 
0 < u < oo,  0 < w < u,  the  electron- electron  contribution  vanishes. 


f -I! 


dudww^^f  (u)  f (w)  = 0 (66) 

o o 


From  the  definition  of  inelastic  collision  rate  given  by  (49),  the  power 
balance  equation  (60)  can  be  written,  finally,  as 

+ e ^U+)  = eu  dng/dt  + 

oo 

ene(2e/m)^2  J*  du  u2«^(2m/M)N f^  + (kT/e)  dfQ/du  J 

o 

+ neS(eu«)  [ Ncr  < vQa(u))-N«<vQ.a(u»] 
a 

co 

+ n e(2e/m)*^  ^ N f du  u^  f (u)  Q (u)  (67) 

e *—i  a J o a 

a o 

where  the  sum  over  "a"  is  over  all  attachment  and  recombination  pro- 
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cesses,  which  are  to  be  implicitly  excluded  from  the  sum  over  all  the 
other  inelastic  collisions  a.  is  the  total  rate  of  creation  of  electrons 
with  positive  energies, 

oo 

Sb  = J du  S (u)  (68) 

o 

and  an  average  deposition  energy  <U+>  has  been  defined  as  a weighted 
average  of  the  external  ionization  source  function  S(u)  as 


oo 


oo 


<U+>  = j*  du  u S(u)  / J"  du  S(u). 


(69) 


The  typical  source  of  external  ionization  is  a high  energy  electron  beam. 
The  mathematical  description  of  the  energy  deposition  in  the  gas  is,  in 
general,  quite  complicated  and  several  investigators  have  developed  ex- 
tensive computer  programs  for  such  analysis.  It  is  beyond  the  scope  of 
the  present  work  to  describe  the  physics  of  e -beam  ionization,  so  it  is 
not  possible  to  indicate  how  the  form  of  the  excitation  source  function 
S(u)  is  determined.  If  the  origin  of  the  source  ionization  is  an  external 

e -beam  of  current  density  J (A/cm^),  then  the  function  S(u)  could  be 

b + 

expressed  in  terms  of  cross  sections  <7\  (u)  as 


S(u)  = (Jb/e)£<r+(u)N., 


(70) 


where  the  sum  is  over  all  species  i ionized  by  the  beam.  The  total  de- 
position of  e -beam  energy  associated  with  ionization  is  then 


= <U+>][V  N.Jb  .JbdV+/dx, 


<?D 
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where  <r  + = J du0\+(u),  and  an  effective  voltage  drop  has  been  defined 
as  dV+/dx  = <U+>/£  in  terms  of  a "mean  free  path"  JL  for  ionization, 


t 


1 


(72) 


A similar  term  dV  /dx  occurs  for  e -beam  deposition  corresponding 
to  excited  states. 


Numerical  Approach 

The  numerical  solution  of  the  Boltzmann  equation  (41)  is  carried  out  by- 

reducing  the  differential  equation  to  a finite  difference  equation,  defined 

over  a uniform  grid  of  equally  spaced  energy  values  (0  = u^,  u^,  u^,  . . . 

. ,.u.  = u ),  where  the  maximum  energy  u satisfies,  typically, 

M max  max 

u > 5u.  This  results  in  a large  set  of  coupled  algebraic  equations 
max 

(typically  M > 500)  which  could  be  solved  (at  least  formally)  by  a matrix 
inversion.  However,  a straightforward  matrix  inversion  is  not  practi- 
cal computationally  or  from  the  standpoint  of  computer  core  storage  that 
would  be  required.  For  example,  a typical  energy  range  of  (0,  20)  eV, 
subdivided  into  a mesh  of  M = 500  to  give  an  energy  resolution  of  0.04 
eV,  would  require  storage  and  manipulation  of  500  x 500  matrices,  each 
of  which  would  contain  250,  000  ~ 1,  000,  000g  elements.  To  reduce  the 
problem  of  core  storage  and  matrix  inversion,  an  iterative  approach 
that  utilizes  the  advantages  of  sparse  matrices  (i.  e.  matrices  whose 
elements  are  mostly  zeros)  has  been  developed.  The  essence  of  the 
present  technique^  is  the  reduction  of  the  system  of  linear  difference 
equations  to  a form  which  involves  the  inversion  of  only  tridiagonal  ma- 
trices, for  which  computation  time  and  core  storage  are  substantially 
reduced. 


With  some  recent  exceptions,  most  of  the  work  which  has  been  done  to 


implement  the  numerical  solution  of  the  Boltzmann  equation  has  been 
associated  with  problems  for  which  the  excited  state  population  densi- 
ties are  negligible,  and  for  which  the  super  elastic  collision  terms  can 
therefore  be  omitted.  In  that  case,  electrons  always  lose  energy  in 
collisions,  and  for  any  given  energy  u the  inelastic  collision  term  on 
the  right  hand  side  of  (41)  involves  only  fQ(u)  and  fQ(u  + ua)--i.  e.  the 
values  for  the  distribution  function  f evaluated  at  energies  equal  to 
and  greater  than  u.  Therefore,  one  of  the  most  common  approaches 
to  the  numerical  solution  of  the  Boltzmann  equation  has  been  to  assume 

a fixed  value  for  f (u  ),  integrate  (41)  backward  from  u = u to 
o max  max 

u = 0,  and  finally  to  impose  the  normalization  condition  (48).  If  the 
superelastic  terms  are  retained,  then  terms  containing  fQ(u  - ua)  oc- 
cur, and  the  backward  integration  scheme  fails.  Attempts  to  extend 
this  technique  to  the  superelastic  case  have  generally  encountered  dif- 
ficulties and  limited  success. 

jjc  -■  5 

Under  highly  excited  conditions  typical  of  laser  plasmas  (N  /N  > 10  ), 

superelastic  collision  terms  must  be  retained.  Their  inclusion  presents 
no  difficulties  for  the  present  technique,  which  is  based  upon  an  itera- 
tive solution  of  the  set  of  linear  difference  equations  obtained  from  (41) 
and  reduced  to  tridiagonal  form.  The  algorithm  developed  is  extremely 
fast,  and  execution  speed  is  optimized  when  a good  initial  guess  for  f 
is  available.  Thus,  the  present  approach  is  especially  suitable  for  a 
coupled  analysis  in  which  the  excited  state  population  densities  and  dis- 
charge voltage  change  as  a function  of  time,  since  the  previous  elec- 
tron distribution  (at  some  time  t)  provides  a good  initial  guess  for  ob- 
taining the  updated  distribution  (at  time  t + At).  As  an  example:  for  an 
electron  grid  of  500  points,  with  15  inelastic  collision  processes  and 

electron- electron  scattering,  the  electron  distribution  function  f con- 

_ 3 ° 

verged  to  an  accuracy  < 2 x 10  after  ~8  iterations  in  ~*>0.7  CF  sec- 


onds  on  the  CDC  6600  computer  (from  an  initial  guess  for  f (u)  given 

- eu/kT  ° 

by  e*  , for  example).  The  numerical  techniques  for  the  solution 

of  Eq.  (41)  will  be  described  in  detail  below. 

The  functional  values  of  the  distribution  f (u)  defined  over  the  uniform 

o 

grid  of  energy  values  (0  = u^  u2*  u3»  • • • » ^ = um  ) shown  in  Fig. 
2.7  will  be  denoted  by  (f^,  f^,  iy  ....  f^),  and  the  spacing  of  the  en- 
ergy grid  will  be  denoted  by  Au: 


M = MESH  + 1 

Au  = u /MESH 
max 

(73) 

u^  = (i  - 1)  Au 


where  i = 1,  2,  3,  ....  M. 


Fig.  2.7  : Uniform  energy  grid. 


If  Eq.  (41)  is  integrated  from  (u.  - Au/2)  to  (u^  + Au/2)  for  i = 2,  3,  4, 
M,  with  the  tacit  boundary  condition  that  f^+j  = 0,  the  result  is 
a system  of  (M-  1)  coupled  algebraic  equations: 
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All  of  the  integrals  which  occur  in  Eq.  (74)  will  be  approximated  by 
the  values  of  the  integrand  taken  at  the  midpoint  u^.  The  values  of  fQ 
and  dfQ/du  which  occur  in  Eq.  (74)  at  the  midpoint  energies, 

u.*  = u.  ± Au/2  (75) 

ii 

are  approximated  by 

fQ(u.  * Au/2)  = fju.*)  = (f.±1  + f.)/2  (76.1) 

f '(u.  ± Au/2)  = f '(u.*)  = - f . ) / A u . (76.2) 

O 1 O 1 lx  1 1 

For  convenience,  we  define  the  following  functions: 

p(u)  = (E2/3)u/<NQ  (u)>  + (kT/e)  u2<(2m/M)NQ  (u)  > 

m m 

+ (2r/3)(q4/e2)nelnAP(u)  (77) 

q(u)  = u2<(2m/M)NQ  (u)  > + ( 2 7r/3)  (q4/e2)  n InA  Q(u)  (78) 

m e 
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in  terms  of  which  the  total  "current  density"  can  be  expressed  as 


J(u)  = J (u)  + J (u)  + J (u) 
i el  ee 

_ _ /_.vi/2r  , „ df« 


= - ne(2e/m)1/2^p(u)  + q(u)  fQ(u)  J (79) 

With  the  above  definitions  and  approximations,  Eq.  (74)  becomes 

‘ p(ui+)(fi+i  * fi,/Au  + P(ui")(fi  ’ fi-i,/Au  ’ ^ui+H£i+i  + V/2 

+ q(u')  (f  + f )/2  + (2e/m)‘1/2Au  . -fi.  - 5.  Av.  +j!°)/Au 

111-1  Lne  1 1 ne  l»z  1 nft 

1 i“E[<“i  + tta)£o(ui+uo)NaQa<"i  + “a>  - ViN«  .> 

a 

+ <Vu0)fo<Vua)N*Q_a(Vua)-  “i£iNaQ.0<“i>]  <S0» 

where  = S(u.).  For  convenience,  define  the  quantities 


Pi  = P*Ui  ^ = P^Ui 


q.*  = qfu.*)  = q(u.  ±Au/2) 


The  5-function  in  Eq.  (41)  and  (74)  is  approximated  by  a rectangular 
spike  of  thickness  Au  and  amplitude  1/Au.  Separate  the  collision  term 
in  Eq.  (80)  into  "diagonal"  (i.e.,  containing  fQ(u.)  * f ) and  "off-dia- 
gonal" elements  (i.e.  containing  f (u.  ± ua)  ) and  use  the  detail  balance 
relations  (33)  to  define 

Di  = £[uiNaQa<ui>  + (ui  + ua)NaQa(ui  + ua)]  (82> 
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and 

Ti  = 23  [(ui  + ua)f0<ui  + ua,NaQa(ui  + ua) 

a 

+ Uifo(ui  ‘ u«)Na  Q«(ui>  ] <83) 

Eq.  (80)  can  be  rearranged  to  give 


a.  . , £ . , + a.  . £.  + a.  ...  f.  , = b., 

i,i-l  i-l  i,ii  i,i+l  l+l.  i 


(84) 


i = 2,  3,  ...»  M,  to  form  a system  of  (M-  1)  equations  for  M unknowns 
(f^,  iy  iy  . ..,  fj^),  with  coefficients  defined  by 


a.  . , = p.  Mu  - q.  /2 

l,  l-  1 i l 


a 


i,  i+1 


p*Mu  + q*/2 


a. ^.  = - (p.+  + p 7 ) / A u + (q.+  - q.")/2 

- Au  [ u.1/2  (2e/m)  ^ n^/n^  + ^i  j (85) 

b.  = - Au  [Ti  + (2e/m)_1/2  [S./ne 

+ Si.2,1'i  + So/n«)/Aul]  ’ 


n 


for  i = 2,  3, . . . , M.  Note  that  the  diagonal  terms  have  been  placed  on  the 
LHS,  and  the  off-diagonal  terms  on  the  RHS  of  Eq.  (84).  There  is,  obvious- 
ly, considerable  latitude  in  the  definition  of  coefficients  (e.g.,  Eq.  (85))  for 
the  reduction  of  the  Boltzmann  equation  to  the  form  (84).  The  essential  cri- 
terion is  that  any  such  choice  be  suitable  for  numerical  implementation  by 
a convergent  and  computationally  efficient  algorithm  for  solution. 
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The  set  of  (M-  1)  equations  (84)  can  be  completed  with  an  Mth  equa- 
tion in  a variety  of  ways.  The  simplest  approximation,  which  has  been 
found  to  be  quite  satisfactory,  is  to  take  ^ After  the  iterative 

calculation  of  f converges,  the  magnitude  is  adjusted  to  satisfy  the  re- 
quired normalization  condition  (48).  Thus,  Eq.  (85)  is  supplemented  by 


a 


1.  1 


1.0 


(86) 


and  the  complete  set  of  M equations  in  M unknowns  can  be  written 


b3  (87) 


aM,  M-lfM-l  + aM,  M fM 

in  which  form  they  resemble  a linear  tridiagonal  system.  It  should  be 
kept  in  mind,  however,  that  they  are  neither  linear  or  tridiagonal;  the 
electron- electron  scattering  contribution  makes  p(u)  and  q(u)  (and  there- 
fore, the  coefficients  a.  a.,  i±j)  dependent  upon  fQ(u),  and  the  vector 
bj  contains  off-diagonal  elements  fQ(u.  ± \ia).  Formally,  the  solution 
of  the  finite  difference  equations  could  be  obtained  (in  the  absence  of 
electron- electron  collisions,  anyway)  by  direct  matrix  inversion,  al- 
this  is  not  practical  computationally  or  from  the  standpoint  of  computer 
storage  requirements.  Instead,  the  set  of  equations  has  been  explicitly 
separated  into  the  form  of  Eq.  (87)  in  order  to  implement  an  iterative 


al,  lfl 


a2,  lfl  + a2,  2f2  + a2,  3f3 


a3,  2f2  + a3,  3f3  + a3,  4f4 
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technique  which  exploits  the  tridiagonal  structure.  Eq.  (87)  can  be 


written  in  matrix  form  as 


where 


A(f)f  = b(f) 


A = 


*1.1  *1.2 

2,  1 2|  2 2,  3 

a a a 
3,  2 3,  3 3,4 


j M-  1,  M 

L aM,  M-  1 aM,  M 

and  f = (fj,  fy  iy  ....  fM)T,  b = (b^  b^,  b^,  bM)T.  The  expli- 
cit dependence  upon  £ of  the  matrix  A(f)  and  the  vector  b ( £ ) has  been 
expressed  in  Eq.  (77),  (78),  (81),  and  (85).  To  solve  the  (nonlinear) 
equation  (88),  consider  a sequence  of  vectors 


{ f 1 } = f*  ° 7n  fn+1 

v * J * | 1 , • • • , A , A ,••• 


rn  -rn+1 


defined  by  the  iterative  equation 


A(fn)fn+1  = b(fn) 


fn+1  = a(  f n)  _1b(  rn), 


with  the  initial  vector  f chosen  arbitrarily.  If  the  sequence  of  vec- 
tors { f defined  by  the  recursion  relation  (91)  converges,  it  must  ne- 
cessarily converge  to  the  solution  of  Eq.  (88),  independent  of  the  initial 
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guess  chosen  for  f . Since  A is  tridiagonal,  it  can  be  numerically  in- 
verted very  rapidly  by  using  elementary  row  manipulations  downward 
(upward)  to  annihilate  the  lower  (upper)  diagonal,  followed  by  back  sub- 
stitution upward  (downward).  Furthermore,  the  computer  core  require- 
ments are  substantially  reduced,  since  only  three  diagonals  (rather  than 
a full  MxM  array)  of  elements  need  be  stored  and  manipulated. 

To  implement  the  iterative  technique,  an  initial  guess  for  f ° is  chosen 
and  a sequence  { f } of  successively  refined  approximations  is  genera- 
ted using  Eq.  (91).  An  acceptable  approximation  f n+*  to  the  solution 
of  Eq.  (88)  is  assumed  to  have  been  attained  when  the  following  conver- 
gence criterion  is  satisfied: 

Max.  | (f“+1-f“)/f"  <e  (92) 


The  solution  f to  Eq.  (88)  is  typically  exponential,  varying  over  sever- 
al orders  of  magnitude  in  the  range  (O.u^^).  Therefore,  Eq.  (92)  has 
been  weighted  in  such  a way  that  the  maximum  relative  (rather  than  ab- 
solute) change  in  ail  of  the  components  T of  the  solution  between  two 
successive  iterations  must  be  less  than  € for  convergence  to  be  de- 
clared. Eq.  (92)  is  not  rigorous,  and  probably  requires  mathematical 
assumptions  about  uniform  convergence;  however,  it  is  physically  rea- 
sonable and  numerically  convenient.  (To  illustrate  the  need  for  caution 
in  applying  Eq.  (92),  consider  the  sequence  {S^}  of  partial  sums,  = 

1 + (1/2)  + (1/3)  + ...  + (1/n).  For  any  fixed  e and  for  n sufficiently 
large,  j)/Sr  = l/nSn<  e.  Although  the  condition  (92)  predicts 

convergence  at  some  point,  the  sequence  (Sn)  is  known  to  diverge.) 

Note  that  all  of  the  off-diagonal  elements  of  the  inelastic  collision  term 
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have  been  included  in  the  vector  b,  defined  in  Eq.  (85).  Terms  involv- 
ing an  energy  (u.iuff)  are  evaluated  at  the  nearest  mesh  point,  or  are 
discarded  if  the  energy  lies  outside  of  the  range  [0,0^^]  chosen  for  the 

calculations.  (It  is  important,  therefore,  to  choose  a value  for  u that 

max 

is  suitably  large  to  minimize  computational  inaccuracies  that  may  result 

from  such  numerical  "sinks".  On  the  other  hand,  u cannot  be  chosen 

max 

too  large  without  sacrificing  energy  resolution.  ) 

There  are  (at  least  in  principle)  an  infinite  number  of  ways  in  which  Eq. 
(80)  can  be  represented  by  a recursive  sequence  of  equations  such  as  (91). 
There  is  no  a priori  guarantee,  in  general,  that  an  iterative  procedure 
based  upon  an  arbitrary  definition  of  the  coefficients  for  the  recursive 
scheme  of  Eq.  (91)  will  necessarily  converge.  In  the  course  of  develop- 
ment of  the  present  algorithm,  some  variations  in  the  definition  of  the 
coefficients  a^  . and  b^  were  explored.  For  example,  in  an  attempt  to 
optimize  performance,  the  definition  (85)  was  modified  by  shifting  a frac- 
tional part  £ of  some  of  the  diagonal  terms  in  a.  . to  the  term  b^  on  the 
RHS  of  Eq.  (84).  It  was  found  that,  in  some  cases,  the  rapidity  of  con- 
vergence could  be  slightly  enhanced  by  choosing  £ > 0,  although  problems 
were  often  incurred  for  large  values  of  £ ( > 0.4)  for  which  convergence 
was  sometimes  unsuccessful.  It  has  been  demonstrated  by  extensive  nu- 
merical tests  that  the  system  of  equations  defined  by  Eq.  (85)- (91)  do  con- 
verge (usually),  and  that  they  are  quite  efficient  computationally.  For  the 
present  analysis,  therefore,  the  definition  (85)  has  been  adopted  for  the 
Boltzmann  subroutine. 

The  coefficient  a^  . defined  by  Eq.  (85)  contains  the  parameter  ne,  whose 
value  is  not  known  a priori.  The  first  integral  of  the  Boltzmann  equation, 
Eq.  (56),  expresses  ne  in  terms  of  the  secondary  ionization,  electron  at- 
tachment and  recombination,  and  external  source  creation  rates.  The 
frequencies  V ^ and  cannot  be  calculated  without  knowledge  of  the  elec- 
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ledge  of  the  distribution  function  fQ(u).  Therefore,  a self-consistent 
iterative  approach  was  taken  to  solve  this  problem. 

An  initial  estimate  (which  may  be  zero)  is  assumed  for  rie,  and  the 
Boltzmann  equation  is  solved  recursively  (as  described  above)  until 
the  convergence  criterion  of  Eq.  (92)  is  satisfied  (with  an  initially 
coarse  choice  of  e,  e.g.,  €~.01).  The  resulting  distribution  fQ(vi)  is 
then  used  to  calculate  the  secondary  ionization,  attachment,  and  re- 
combination rates  which  can  be  used,  along  with  the  external  creation 
rates,  to  obtain  an  improved  estimate  of  rie  from  Eq.  (56).  The  pro- 
cedure is  repeated  until  the  value  for  ng  converges  to  some  accept- 
able accuracy.  After  convergence  relative  to  the  coarse  parameter 
of  accuracy  is  obtained,  further  refinement  can  be  obtained  by  making 
the  parameter  6 smaller. 

The  implementation  of  these  numerical  algorithms  is  accomplished 
by  a FORTRAN  subroutine  BOLTZ,  to  be  described  in  more  detail  in 
Sec.  3.3  below.  In  addition,  the  subroutine  is  extensively  document- 
ed with  internal  COMMENT  cards  (cf.  listing  in  Vol.  II).  The  pres- 
ent version  of  Program  LASER  uses  the  Boltzmann  analysis  only  for 
the  situation  of  an  electric  discharge,  although  in  principle,  the  for- 
mulation should  also  be  applicable  to  the  case  of  direct  e"-beam  exci- 
tation only. 


! 
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2.4  External  Driving  Circuit  and  Ionization  Sources 

A typical  excitation  scheme  for  high  power  electrical  lasers  employs 
a capacitive  electric  discharge,  sustained  by  an  external  source  of  io- 
nization such  as  a high  energy  electron  beam.  In  some  cases,  the  ad- 
ditional enhancement  of  electrical  pump  power  provided  by  a discharge 
is  not  justified  in  terms  of  the  additional  complexity  or  problems  that 
result,  and  direct  excitation  by  e -beam  alone  may  be  used.  The  pre- 
sent analysis  includes  external  sources  of  ionization,  as  discussed  in 
Sec.  2.3  for  the  electron  kinetics,  and  couples  the  kinetic  equations  to 
the  voltage  and  current  equations  for  an  external  driving  circuit. 

The  model  for  the  external  driving  circuit  is  shown  in  Fig.  2.8.  A 
plane  parallel  discharge  of  area  A,  with  anode  and  cathode  separated 
by  a distance  d,  is  driven  by  an  external  RLC  circuit  with  a capacitor 
initially  charged  to  a high  voltage.  This  discharge  is  sustained  by  an 
external  e -beam  whose  temporal  current  density  J^( t)  can  be  speci- 
fied as  an  arbitrary  function  of  time.  Initially,  there  is  no  voltage  a- 
cross  the  discharge  because  of  the  external  circuit  inductance.  As 
the  secondary  electron  density  in  the  gas  changes,  the  plasma  conduc- 
tivity (and  hence  the  discharge  impedance)  will  change.  Thus,  the 
plasma  kinetics  analysis  is  directly  coupled  to  the  circuit  equations 
which  define  the  instantaneous  voltage  V^(t)  (and  thus,  E(t)/N)  across 
the  discharge.  Although  the  present  model  for  the  external  driving 
circuit  is  the  simplest  possible,  it  neverthless  provides  for  the  de- 
scription of  basic  but  important  phenomena  associated  with  pulsed  ex- 
citation of  the  medium.  If  necessary,  more  sophisticated  circuits 
(e.g.,  pulse- forming  networks  to  better  match  impedance^)  could  be 
incorporated  into  the  analysis. 
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EXTERNAL  DRIVING  CIRCUIT 


Plasma  conductivity*,  a • ne(t)ep(t) 
Plasma  impedance:  Rd(t)  - d/Aa(t) 

lIL  + (r  + Rj(t))  i + 0 

dt  d C 


Fig.  2.8:  External  RLC  driving  circuit  excites  gas  in  a 

plane- parallel  discharge  of  area  A,  separation  d 


6 


The  circuit  equations  will  be  formulated  as  first  order  differential 
equations  for  the  charge  q(t)  and  current  i(t)  in  order  to  incorporate 
them  on  an  equal  basis  with  the  molecular  kinetic  equations  (4),  the 
full  set  of  which  are  integrated  numerically  using  the  Gear^  technique. 
Thus, 


dq/dt  = i(t)  (93.1) 

Ldi/dt  = -q/C  - (r  + Rd(t))i  (93.2) 

where  L = inductance,  r = external  resistance,  and  C = capacitance. 
If  the  case  L = 0 is  desired,  only  one  circuit  equation  is  required: 

dq/dt  = q/(r  + Rd(t))C  (94) 

The  discharge  impedance  Rd(t)  is  a function  of  time,  given  by 

Rd(t)  = d/A<r(t),  (95) 

where 

o'(t)  = ne(t)e/i(t)  (96) 

is  the  plasma  conductivity  as  a function  of  the  instantaneous  electron 
density  ng(t)  and  the  mobility /u(t).  The  voltage  across  the  discharge, 


Vd(t)  = i(t)Rd(t) 


(97) 


determines  the  value  of 


E(t)/N  = Vd(t)/dN 


(98) 


If  Eq.  (93.2)  is  multiplied  by  i(t),  an  expression  for  the  conservation 
of  energy  (power)  is  obtained: 


I 


d/dt  (Li2/2)  + i2r  + (Ad)o-E2  = - d/dt  (q2/2C). 


(99) 


The  physical  interpretation  of  Eq.  (99)  is  that  the  rate  of  change  of 
stored  inductive  energy  (Li2/2),  plus  power  dissipated  in  the  external 
resistor  (i  r)  and  deposited  into  the  discharge,  is  equal  to  the  rate  of 
change  of  stored  energy  in  the  capacitor  (q2/2C).  The  discharge  power 
density 


o-E2  = ng(t)e^(t)E2  (100) 

has  occurred  previously  in  Eq.  (67)  for  electrical  power  partitioning. 
In  the  following  section,  the  overall  power  balance  equation  for  the 
complete  coupled  system  will  be  discussed. 

2.  5 Conservation  of  Energy 

In  the  foregoing  sections,  2.  1 - 2.4,  equations  for  power  balance  have 
been  developed  for  each  component  of  the  coupled  analysis.  It  is  now 
appropriate  to  combine  these  results  and  to  formulate  a power  balance 
relation  for  the  complete  system.  The  physical  interpretation  of  the 
parameter  AE(ot)  for  each  reaction  determines  its  physical  signifi- 
cance in  the  power  balance  equation.  To  begin,  therefore,  some  fur- 
ther discussion  devoted  to  interpretation  of  the  terms  of  Eq.  (5)  for 
specific  types  of  kinetic  reactions  is  required.  This  will  provide  the 
opportunity  to  describe  the  typical  types  of  reactions  permitted,  and 
various  conventions  or  restrictions  on  syntax. 
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Electron-  Beam  Collisions 

It  is  beyond  the  scope  of  the  present  analysis  to  give  a rigorous  treat- 
ment of  electron-bei..—  energy  deposition  in  a gas,  although  provisions 
have  been  made  for  the  inclusion  of  high-energy  electron  collision  pro- 
cesses in  the  molecular  kinetic  reaction  scheme  in  terms  of  a simple 
rate  constant  formulation.  Excitation  and  ionization  of  the  molecular 
species  X by  collision  with  high-energy  electrons  T , 


$ 

X + e"  — X + e'  + AE 


(101.  1) 


X + e"  — X+  + e"  + e"(u)  + AE  (101.2) 


* + 

are  described  in  terms  of  cross  sections  <r  and  ar  (u),  respectively. 

It  is  assumed  that  the  medium  is  isotropically  excited  and  "thin",  so 
that  the  approximation  of  single  scattering  may  be  made.  Adjustments 
for  the  effects  of  multiple  scattering  can  be  made  in  a qualitative  way 
by  multiplying  the  cross  sections  by  an  effective  scale  factor  if  some 
a priori  knowledge  of  electron-beam  deposition  as  a function  of  electron 
energy,  gas  mixture,  and  gas  pressure  is  available.  In  terms  of  the 
e’-beam  current  density  J (A/cm  ),  the  rates  for  reactions  (101)  are 


-d[x]/dt  = d[x*]/dt  = (Jbtr*/e)[x]  (102.1) 

-d[x]/dt  = d[X+]/dt  = (Jbcr+/e)[x]  (102.2) 

where  or+  = f duo^u).  The  energies  AE,  as  defined  by  Eq.  (6),  for 

J * + 
the  two  reactions  (101)  are  AE  = -eu  and  AE  = -eu  respectively, 

where  u*  and  u+  are  excitation  and  ionization  energies  (eV).  Thus, 

the  contribution  of  e"-beam  collisional  excitation  to  the  molecular  power 
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balance  Eq.  (5)  becomes 


I 


( 

i 


£(.u.*H[x’]/dt  --  JbJ.V*[X.] 
i i 

(103.  1) 

2<.„.+)d[x+]/d.  = dbX<  <[x.] 

i i 

(103.2) 

where  the  sums  are  over  all  species  which  are  excited  or  ionized  by 
the  e -beam.  Strictly  speaking,  the  reverse  process  for  reaction 
(101.  1)  should  also  be  included,  although  its  contribution  can  usually 
be  neglected  in  Eq.  (103.  1).  When  Eq.  (67)  for  electron  power  balance 
is  added  to  Eq.  (5)  for  the  molecular  power  balance,  it  is  apparent 
that  the  contribution  of  the  e -beam  will  be  given  by 

co 

JbdV/dx  = <r*  + u.+  <r.+  + j*  duu<r*(u)][X.] 

i o 

* \ 2[ai*  \ + | \ + <u+>j«r+  ] [x.].  (104) 

i 

where  dV/dx  is  an  effective  voltage  drop  for  the  beam,  and  the  aver- 
age energy  <U+>  deposited  into  the  electron  gas  (defined  originally  in 
Eq.  (69))  is  given  by  the  weighted  average, 


oo 

<U+>  = £[X.]J  duu<r+(u).  (105) 

i i o 

Cross  section  data  for  high-energy  electron  beam  excitation  has  been 

23 

given  by  Berger  and  Seltzer  . If  multiple  scattering  effects  are  im- 
portant for  the  physical  problem  to  be  analyzed,  the  Berger- Seltzer 
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cross  sections  can  be  appropriately  adjusted  by  multiplying  by  some 

scale  factor  to  produce  an  effectively  larger  deposition.  The  present 

analysis  requires,  as  input  rate  constants,  the  total  scattering  cross 
* + 

sections  o\  and  cr.  for  the  e -beam  collision  processes.  In  the 
i i * 

case  of  excited  state  production,  the  information  provided  by  °\  is 
sufficient.  However,  e -beam  ionization  processes  must  be  included 
in  the  source  term  (70)  to  the  Boltzmann  equation,  and  therefore,  the 
detailed  energy  dependence  of  the  cross  section  o\+( u)  is  required. 
For  simplicity,  it  shall  be  assumed  that  the  shape  of  the  energy  depen 
dence  is  the  same  for  all  processes,  so  that  the  source  term  given  by 
Eq.  (70)  can  be  written 


where  s(u)  is  a normalized  shape  function  satisfying 


Thus,  the  average  energy  <U  > deposited  as  electron  kinetic  energy 


The  total  rate  of  creation  of  secondary  electrons  by  e -beam  ioniza- 

.3 

tion  (cm  /s)  is  thus 


— 


Sb=JduS(u)  = (Jb/e)^o-.+  N.  Jdus(u)  = (Jb/e)  ^<r.+  N.,  (109) 


as  given  by  Eq.  (68),  and  the  energy  deposition  rate  (i.e.,  power  den- 
3 

sity,  W/cm  ) is 


j*du(eu)S(u)  = eSQJ*duus(u)  = e^U^S^, 


(110) 


where  e = 1.602  x 10  7 J/eV.  For  convenience,  the  normalized  shape 

function  s(u)  is  defined  (by  input)  to  be  a square  wave, 


u < u , 

cL 


s(u)  = < l/fi^-u^,  ua^u^ub, 


(111) 


and  thus, 


<U+>  = (u^  + *1^/2. 


(112) 


Because  the  present  approach  to  e -beam  deposition  is  phenomenologi- 
cal with  no  a priori  knowledge  of  energy  dependence,  the  simplified  as- 
sumptions for  s(u)  must  suffice.  Furthermore,  it  should  also  be  noted 
that  the  present  formulation  of  e -beam  excitation  and  ionization  colli- 
sions does  not  provide  an  explicit  mechanism  for  heating  the  gas  by  di- 
rect deposition  into  molecular  kinetic  energy.  The  phenomenological 
description  of  e -beam  deposition  should  be  sufficiently  flexible  to  be  ap- 
plicable to  typical  situations  of  interest.  Where  the  approximations  are 
inadequate,  the  analysis  must  necessarily  be  extended  as  required. 
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In  addition  to  e -beam  ionization,  there  are  a variety  of  kinetic  mecha- 
nisms which  involve  the  creation  of  secondary  electrons.  It  is  appro- 
priate to  make  some  comments  here  with  regard  to  the  manner  in  which 
such  processes  are  to  be  included  in  the  Boltzmann  analysis,  and  to  dis- 
cuss the  resulting  consequences  for  the  power  balance  equation. 

In  general,  any  kinetic  reaction  which  involves  the  creation  of  secon- 
dary electrons  must  be  accounted  for  in  the  source  term  Sext(u)  for  the 
Boltzmann  equation.  The  contribution  of  e -beam  ionization  to  S(u)  is 
given  by  Eq.  (106).  Other  electron  creation  processes  could  include, 
for  example, 

Photoionization:  A + I(v)  — A+  + e 

Photodetachment:  A + I(V)  -•  A + e 

* + - 

Penning  Ionization:  A + B A + B + e 

Since  there  is  no  a priori  knowledge  of  the  kinetic  energy  distribution 
of  the  created  secondary  electrons,  it  is  convenient  to  assume  that  they 
are  all  created  at  zero  energy  and  that  the  remaining  energy  AE  is  con- 
verted to  molecular  kinetic  energy.  (This  is  analogous  to  the  treatment 
of  electrons  created  by  secondary  ionization.)  The  contribution  from 
such  processes  to  the  source  term  is  proportional  to  a 8 -function: 

Sext  (u)  = (Jb/e)  ^°’i+Nis^u)  + S (“)  [Iaph  [a]  /hV  + 

i 

+ I<rpd[A']/hV  + kpl[A*][B]+  ...]  (113) 

The  8 -function  term  makes  no  contribution  to  electrical  power  balance. 
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Secondary  Electron  Collisions 


For  a binary  electron-molecule  collision, 


X+e  — X + e + AE 


the  energy  AE  is  given  by  AE  = - eu  , where  u is  the  inelastic  energy 
loss  (eV)  for  the  process.  Thus,,  the  contribution  of  secondary  electron 
collisions  to  the  right  hand  side  of  the  molecular  power  balance  Eq.  (5) 


£r(«)AE(«)  = ne^(eua)  [Na<VQa>  - Na<vQ_a>]  (114) 


Note  that  the  sum  on  the  right  hand  side  of  Eq.  (114)  is  not  restricted  to 
exclude  attachment  (or  recombination)  processes,  as  was  the  case  in  Eq. 
(67)  for  electrical  power  balance.  For  attachment  and  recombination, 
denoted  by  "a",  the  difference 


N e(2e/m)^  f duu^f  (u)  Q (u)  - n N (eu  )^vQ  > 
a J oa  e a a a 

o 

oo 

= n N e(2e/m)*^  f duu(u-u  )f  (u)Q  (u)  (115) 

e a J a o a 


can  be  interpreted  as  (a  contribution  to)  the  kinetic  heating  of  the  molecu- 
lar gas  mixture.  For  example,  in  a "two-body"  recombination  process 
such  as 


X+e  - X 


(116) 
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to  kinetic  energy  of  X and  to  recombination  radiation  energy  - hV^ 
of  a photon  (which  has  been  tacitly  suppressed  in  (116) ).  Since  u^  ( = 

- AE)  = - u+  for  reaction  (116),  the  resulting  power  density  is 

I 

oo 

P = eng[X+]  (Ze/m)1^  J du  u ^(u  + u+  - uy  ) + uyJ  fQ(u)  Qr(u),  (117) 

o 

where  Q (u)  is  the  recombination  cross  section.  The  first  term  repre- 
r 

sents  conversion  of  energy  to  molecular  kinetic  energy,  and  the  second 
to  radiation.  In  the  subsequent  discussion,  we  shall  neglect  the  distinc- 
tion, and  simply  interpret  Eq.  (117)  as  molecular  kinetic  heating. 


Radiative  Processes 


Both  spontaneous  as  well  as  stimulated  emission  and  absorption  proces- 
ses can  be  included  in  the  present  analysis.  In  addition  to  stimulated 
processes  involving  the  laser  transition,  absorption  of  laser  radiation  by 
other  species  may  also  occur.  There  are  certain  distinctions  and  con- 
ventions that  must  be  observed  in  the  present  analysis,  and  the  corres- 
ponding interpretations  for  the  power  balance  equation  will  be  discussed. 


First  of  all,  spontaneous  emission  photons  (except  for  those  which  build 
up  the  laser  field)  are  neglected  as  independent  species  in  the  reaction 
syntax,  specified  by 

X*  — X.  Cl  I®) 

The  contribution  of  this  reaction  to  the  right  hand  side  of  Eq.  (5)  is  just 

■ S4E  (119) 

sp 
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where  AE  = eu  is  the  energy  defined  by  Eq.  (6),  and  r is  the  spon- 

sp 


taneous  radiative  lifetime.  For  all  of  the  spontaneous  emission  proces- 
ses from  excited  states,  the  terms  of  the  form  (119)  can  be  interpreted 
as  fluorescence  loss.  The  syntax  of  Eq.  (118),  which  explicitly  omits 
reference  to  the  radiated  photon,  is  recognized  by  the  computer  analysis 
to  represent  spontaneous  emission. 


For  the  photon  field  of  the  laser,  however,  the  syntax  for  defining  spon- 
taneous radiation  is 


X — X + hi/. 


(120) 


where  the  emitted  photon  must  be  explicitly  indicated.  The  reason  for 
this  is  that  the  laser  field  is  defined  in  terms  of  the  intracavity  photon 
number  density,  and  laser  photons  are  treated  as  species  on  an  equal  ba- 
sis with  electrons,  ions,  molecules,  and  all  of  the  other  "particle"  spe- 
cies which  occur  in  the  kinetic  reaction  scheme.  To  insure  a parallel 
treatment,  some  further  observations  and  definitions  must  be  made. 


First  of  all,  since  the  volume  of  the  active  gain  medium  is  less  than  the 


cavity  volume  by  a factor  (lWLc),  it  is  appropriate  to  define 


N = (L  /L  )n  . 
ph  c g'  ph 


(121) 


as  an  effective  photon  number  density  in  the  gain  medium  (as  if  the  radi- 
ation were  considered  to  be  confined  to  the  gain  volume).  This  is  logi- 
cal, since  the  power  balance  Eq.  (5)  is  formulated  in  terms  of  the  power 
density  in  the  gain  medium,  and  all  of  the  other  species  are  described 
by  their  population  densities  in  the  gain  medium.  In  terms  of  the  photon 


I 1 
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density  N the  intracavity  intensity  is 


I = (L  /L  )chvN  , 
g c ph 


(122) 


and  the  radiation  Eq.  (13)  becomes 


dNph/dt  = (n/47r)N  Jr  + {WLJc[a(t}  - a^jN 


g c 


thJ  ph* 


(123) 


a(t),  the  instantaneous  net  gain  coefficient,  was  defined  in  Eq.  (10),  and 
is  the  threshold  loss  coefficient,  given  by 

«th  = O + (l/2)ln(l/R)]  /L  (124) 

© 

where  y is  the  intracavity  loss  per  pass  (i.  e.  from  optical  elements), 
and  R is  the  output  coupling  reflectivity.  The  total  contribution  from  out- 
put coupling  and  the  radiative  processes  of  Eq.  (7),  (3),  and  (120)  to  the 
power  balance  equation  (5)  is 


EAd[A]/dt  + EA*d[A*]/dt  + hydN  h/dt 

A A*  P 

+ E^ d [x*] /dt  + Exd[x]/dt  = -£l  - (fl/47r)j  [X*]  hv/r 

- (I/hV)  J <rA[A](EA-  E^+hu)  - I«th  (125) 

A 
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The  interpretation  of  the  terms  of  Eq.  (125)  is  straightforward.  On 
the  left  hand  side,  there  is  the  rate  of  change  of  stored  energy  density. 
The  first  term  on  the  right  hand  side  represents  the  net  fluorescence 
power  loss  from  reaction  (120),  and  is  consistent  with  the  assumption 
that  a fractional  part  (fl/4-jr)  of  the  spontaneous  emission  photons  are 
not  lost,  but  remain  in  the  cavity  to  build  up  the  optical  field.  Refer- 
ence to  Eq.  (8)  shows  that  the  energy  difference  (E  + hV  - E *)  cor- 

A A 

responds  to  the  net  energy  converted  to  kinetic  energy  of  the  species 
A . Thus,  the  second  term  on  the  right  hand  side  of  Eq.  (125)  can  be 
interpreted  as  a kinetic  heating  power  density.  The  third  term  corre- 
sponds to  optical  loss,  due  to  absorption  by  intracavity  optical  elements 
(y),  and  output  coupling  (R). 

Kinetic  Heating 

With  the  exceptions  noted  above,  the  energy  change  AE  defined  by  Eq. 

(6)  for  the  reaction  (1)  is  the  net  energy  converted  to  molecular  kinetic 
energy,  and  thus,  the  physical  significance  of  the  term 

^R(£)AE(/3) 

0 

on  the  right  hand  side  of  Eq.  (5)  is  that  it  represents  power  density  in- 
to kinetic  heating  of  the  molecular  gas.  The  exceptions  which  occurred 
were  the  result  of  the  tacit  suppression  of  some  species  involved  (e.g., 
the  emitted  photon  in  spontaneous  radiation)  or  of  the  kinetic  energy  de- 
pendence of  some  of  the  interacting  species  (e.g.,  secondary  and  high- 
energy  e -beam  electrons).  In  the  case  of  secondary  electrons,  for 
which  an  integration  over  a continuum  of  energies  is  required,  the  con- 
tribution to  the  power  balance  must  be  derived  from  the  Boltzmann  equa- 
tion. 
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Total  Power  Balance  for  the  Coupled  System 


As  a result  of  the  foregoing  discussion,  which  should  clarify  the  physi- 
cal significance  of  the  terms  of  Eq.  (5),  it  is  now  possible  to  construct 
the  equation  for  conservation  of  power  density  for  the  coupled  system  of 
electrons,  molecules,  radiation,  and  external  driving  circuit.  If  Eq. 

(5)  is  added  to  the  expression  given  by  Eq.  (67)  for  eudne/dt,  and  the 
terms  are  rearranged,  the  following  result  is  obtained: 


1 - {n/4ir)l  hv[  X*]/ 


du  u <(2m/M)NQ  > f + (kT/e)df  /du 


crE  + J dV/dx 


The  first  three  terms  on  the  left  correspond  to  the  rate  of  change  of  the 
energy  density  stored  in  the  molecules,  radiation  field,  and  electrons. 
The  sum  over  i excludes  electrons,  but  photons  (which  have  been  shown 
explicitly  in  the  second  term)  could  be  included  in  that  sum  by  defining 


T 


E.  = hv  for  that  "species".  The  fourth  and  fifth  terms  represent  spon- 
taneous fluorescence  loss,  while  the  sixth  is  the  optical  power  loss  by 
intracavity  absorption  and  output  coupling.  The  remaining  terms  on  the 
left  hand  side  have  the  significance  of  kinetic  heating  of  the  molecular 
gas.  The  sum  over  /3  in  the  seventh  term  includes  all  processes  except 
spontaneous  emission,  secondary  electron  collisions,  and  high-energy 
e -beam  collisions.  All  resonant  energy  transfer  processes  (i.  e.  for 
which  AE  = 0)  contribute  nothing  to  the  sum  over  ft,  nor  does  the  stimu- 
lated emission  reaction,  since  AE  = E * - E - hi/  = 0.  The  eighth 
term  consists  of  the  net  contribution  to  molecular  heating  from  recom- 
bination and  attachment  processes,  assuming  that  recombination  radia- 
tion is  neglected.  The  last  term  on  the  left  hand  side  is  the  power  den- 
sity into  elastic  heating  from  electron- molecule  momentum  transfer  col- 
lisions. 

On  the  right  hand  side,  there  is  the  input  electric  discharge  power  den- 
2 

sity  & E and  the  beam  deposition  J^dV/dx,  where  the  effective  e -beam 
voltage  drop  dV/dx  was  defined  in  Eq.  (104).  The  discharge  power 
density  <rE^  can  be  expressed  in  terms  of  the  parameters  of  the  exter- 
nal driving  circuit  by  Eq.  (99). 

In  its  present  formulation,  the  coupled  analysis  assumes  that  the  mole- 
cular temperature  of  the  gas  remains  constant.  This  restriction  could 
easily  be  removed,  if  necessary,  by  an  additional  equation 

dT^^j/dt  = ^ R(/3)AE(j8)  + misc.  heating  (127) 

P terms 

3 

where  is  the  specific  heat  of  the  gas  at  constant  volume  ( J /cm  K), 
and  the  right  hand  side  is  the  total  power  density  into  kinetic  heating. 


74 


The  numerical  solution  of  the  coupled  system  of  molecular  and  electron 
kinetic  equations  proceeds  as  follows.  The  total  pulse  duration  is  divi- 
ded into  a specified  number  (e.g.,  200)  of  (equal)  subintervals  defined 
by  a time  grid  [0,  tJf  tv  . . . , t^,  . . . ].  At  each  of  the  discrete  points  ^ 
in  time,  the  Boltzmann  equation  is  solved  (as  a function  of  the  instanta- 
neous values  E(t^)/N  and  the  excited  state  population  densities  Na  (t^) ) 
to  update  the  electron  excitation  rates  <vQq>  and  plasma  parameters, 
which  are  then  assumed  to  remain  constant  over  the  following  subinter- 

time.  The  molecular  kinetic  (and  circuit)  equations  are 
integrated  over  this  subinterval,  and  in  practice,  the  integration  step 
size  will  usually  be  much  smaller  than  the  coarse  time  intervals  [t^,  i 
over  which  the  electron  kinetics  calculations  are  performed.  The  Gear** 
technique  automatically  adjusts  the  step  size  to  maintain  accuracy  and 


] 


stability. 


The  assumption  that  the  electron  energy  distribution  function  fQ(u)  is  con- 
stant over  these  discrete  subintervals  will  result  in  a power  balance  dis- 
crepancy in  the  calculations.  Assume  that  f*  (u),  <vQtt)k,  pk,  ... 
represent  the  electron  distribution,  excitation  rates,  mobility,  average 
energy,  etc.  obtained  from  solution  of  the  Boltzmann  equation  at  some 
discrete  point  t = tjj.  If  Eq.  (67)  were  to  be  numerically  evaluated  at  a 
later  time  t by  the  approximations  that  fQ(u)  = f^(u),  <vQa>  = <vQ a>k, 

u = uk*  ^ ~ ^ k»  • • • t^ie  right  hand  side  would  contain  an  additional  term 
6P  representing  a power  discrepancy: 


®P(tk»t)  = ene(t)/ikE(t)^  + e^.U+^S^(t)  - eukdne(t)/dt 

oo 

- «e(t)e(2 e/m)1^2  J* duu2<(2m/M)NQm>£f^(u)  + (kT/e)  df*(u)/du  j 
o 

(continued) 
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- ne(t)  ^ (eua)  [Na(tKvQa)k  - N*  (t)<vQ.a>k] 


oo 


- ng(t)  e(2e/m)^2  ^ Na(t)J*  duu^  f^(u)  Qa(u) . (128) 


Since  fQ  (u)  and  the  corresponding  plasma  parameters  and  excitation 

rates  are  recalculated  at  the  point  t , the  longest  period  of  time  over 

k ^ ^ 

which  the  approximation  f (u)  = f (u)  is  sustained  is  the  length  of  the 

o o 


discrete  interval.  At  = t - t . Therefore,  a useful  measure  of  the 

xC*r  1 xC 

accuracy  of  the  calculation  (at  least  insofar  as  it  is  affected  by  errors 
resulting  from  the  discrete  time  partition)  should  be 


«<V‘)  = pVt)/p«„t(t)- 


(129) 


where 


PtQt(t)  = ene(t)|i(t)E(t)  + Jb(t)dV/dx 


(130) 


is  the  total  electrical  input  power  density.  If  Eq.  (128)  is  evaluated  at 
the  time  t = tk,  a measure  of  power  balance  inaccuracy  resulting  from 
the  numerical  solution  of  the  Boltzmann  equation  is  obtained.  (Since  it 
is  only  the  contribution  e<U't">S^)  to  J^dV/dx  that  occurs  in  the  electron 
kinetics  analysis,  the  fractional  power  discrepancy  reported  from  the 
Boltzmann  calculations  is  normally  defined  in  terms  of  total  electrical 
excitation  power  density. 


Pe  = enepE2  + e<U+>Sb. 


(131) 


The  present  formulation  of  the  coupled  system  of  equations  makes  the 
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tacit  assumption  that  the  charge  q(t),  current  i(t),  and  all  of  the  popula- 
tion densities  Na(t)  are  continuous  from  t^  = t^  - 0 to  t^+  = t + 0: 

i(tk_)  = i(tk+)’ 

q(tk")  = q^),  (133) 

Na(\)  = Na(tk+). 

and  thus,  these  quantities  can  always  be  defined  unambiguously  as  func- 
tions of  the  time  t^.  However,  all  of  the  electron  parameters  and  exci- 
tation rates  are  discontinuous  from  t^  to  t^+,  since  they  are  updated  at 
t = t^  by  a new  solution  of  the  Boltzmann  equation.  Furthermore,  the 
derivatives  di(t)/dt,  dne(t)/dt,  and  dNa(t)/dt  for  all  species  a which  are 
involved  in  secondary  electron  collision  processes  will  also  be  discontin- 
uous. At  t = tk~,  just  prior  to  the  Boltzmann  update,  the  electric  field 
is  obtained  from 

Ek  = E(t^)  = il^R^'j/d  = i(tk) / [Aene(tk)  p (tjJ ) ],  (134) 

and  the  fractional  power  discrepancy  s (t^*  defined  by  Eq.  (128)  and 
(131)  is  based  upon  the  assumption  that  Ek  is  the  electric  field.  However, 
the  plasma  conductivity  cr  = en  p (and  thus,  the  discharge  impedance 
Rd  = d/A 
field  at  t = t^+  is 

Ek  * *<tk)  = i(tk)/[Aane(tk)H(t^)]  = (\/#*k)Ek.  (135) 

-bus,  m addition  to  the  power  discrepancy  <5P(tk»  ^ inherent  in  the  nu 
eolation  of  the  Boltzmann  equation,  there  will  also  be  an  imme- 


a)  are  discontinuous  from  ^ to  , and  therefore,  the  electric 
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diate  discrepancy  associated  with  the  discontinuity  in  the  electric  field: 

ap(V*k>  * 8PW  + t(Ek’2  - Ekl 

■ 5p<tk'  *k>  + Ek  t"‘k/“k)2  - »]•  <136> 


where  = p (t^ ) . There  are  three  fractional  power  discrepancies 
which  are  a useful  measure  of  the  accuracy  of  the  calculations.  The  first 
is  the  inherent  power  accuracy  of  the  numerical  solution  of  the  Boltzmann 
equation, 


fk  - 6P<VV/Pe-  <137> 

and  it  is  provided  as  part  of  the  output  summarizing  the  results  of  the  el- 
ectron kinetics  analysis.  The  second  is 

* *(tk.i-V  * sp,tk.rV/p«ot  <138> 

which  represents  the  cumulative  error  which  has  developed  over  the  en- 
tire preceding  subinterval  of  time.  The  third  is 

fk  ’ f(V‘k'  = ^'V'k'^to.  <139> 

which  represents  the  fractional  additional  error  resulting  from  disconti- 
nuity in  the  electric  field  before  and  after  the  Boltzmann  update.  If 
or  are  too  large,  the  length  At  = t^+1  - t^  of  the  subintervals  should 
be  reduced. 
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3.0  COMPUTER  PROGRAM  DESCRIPTION 


I 


The  generalized  laser  kinetics  code  has  been  written  in  FOR  TRAN  IV, 
and  was  developed  for  use  with  the  Extended  FORTRAN  Compiler  of 
the  CDC  6000  and  CYBER  Series  of  computers  (SCOPE  3.4).  If  it  is 
necessary  to  execute  this  program  on  a different  computer  system, 
some  modifications  of  nonstandard  syntax  or  conventions  may  be  re- 
quired. In  addition,  extensive  usage  of  ENCODE  and  DECODE  state- 
ments (for  output  formatting)  makes  the  explicit  assumption  that  the 
basic  "word"  size  is  60-bit,  10-BCD  character  display.  Modification 
of  such  statements  will  be  necessary  if  the  program  is  to  be  adapted 
to  a computer  with  a different  word  size. 

The  program  structure,  flow  diagrams,  description  of  subroutines, 
control  card  commands  and  requests,  external  file  usage,  input  card 
format  and  default  conditions,  available  output  options,  and  illustra- 
tive input  and  output  from  a sample  case  will  all  be  described.  The 
complete  FORTRAN  program  listing  is  given  in  Vol.  II  of  the  present 
report.  In  addition,  Vol.  II  also  contains  a listing  for  an  electron  ki- 
netics analysis  program,  which  makes  use  of  the  same  subroutines. 

3.1  Program  Structure  and  Flow  Diagram 

The  program  structure  has  been  designed  with  the  objective  of  provi- 
ding maximum  input  and  output  flexibility,  while  reducing  mechanical 
programming  tasks  for  the  user  to  a minimum.  A given  kinetic  reac- 
tion scheme  is  specified  by  input  (with  a completely  free  and  flexible 
syntax)  as  a sequence  of  symbolic  reactions  containing  an  arbitrary 
number  of  processes  and  interacting  species.  These  reactions  are 
translated  into  computer- coded  FORTRAN  equations  to  generate  mo- 
lecular kinetics  subroutines,  which  are  then  combined  with  the  master 
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executive  program  LASER  and  (seventeen)  other  subroutines  to  form 
a completely  self-contained  source  program  for  the  coupled  analysis 
based  upon  the  specified  reaction  scheme  (cf.  Fig.  2.3). 


Fig.  3.1  presents  a more  detailed  flow  diagram  of  the  main  program 
LASER.  Initially,  there  may  or  may  not  be  a synthesized  program 
available,  and  the  control  cards  and  data  input  deck  structure  must 
be  compatible  so  that  LASER  can  execute  either  situation.  Initially, 
the  program  makes  an  attempt  to  read  "TITLE"  on  an  attached  file 
TAPE4,  which  would  normally  contain  the  "name"  and  miscellaneous 
rate  constant  data  associated  with  a synthesized  reaction  scheme  of 
a previous  run.  If  an  EOF  is  encountered  on  the  attempt  to  read  file 
TAPE4,  the  program  assumes  that  user  intends  to  generate  subrou- 
tines and  a data  file  for  a new  reaction  scheme,  to  be  specified  by  in- 
put DATA  BLOCK  1.  In  that  case,  the  input  of  DATA  BLOCK  1 is 
processed  by  Subroutine  SYNTH,  which  translates  the  reaction  struc- 
ture into  three  subroutines  (DNDT,  JACOB,  and  LEVELS)  which  con- 
sist of  80-BCD  character  records  on  file  MTAPE  (=  3).  Furthermore, 
as  the  reaction  scheme  of  DATA  BLOCK  1 is  processed,  miscellane- 
ous data  is  written  onto  file  TAPE4.  Diagnostics  are  generated  for  a 
variety  of  error  conditions  which  may  be  encountered.  Although  syn- 
thesis proceeds  by  ignoring  unacceptable  syntax  or  other  errors,  pro- 
tection against  subsequent  execution  of  a faulty  reaction  scheme  can 
be  obtained  by  declaring  such  conditions  to  be  fatal.  After  files  TAPE4 
and  TAPE3  are  constructed,  EXIT  occurs,  the  synthesized  subrou- 
tines are  compiled  and  combined  with  the  other  subroutines,  and  LA- 
SER is  again  executed.  This  time,  TAPE4  contains  reaction  scheme 
data,  and  the  attempt  to  read  "TITLE,  " the  "name"  of  the  synthesized 
program,  is  successful.  Control  transfers  to  Read  DATA  BLOCK  2, 
which  contains  input  electron  cross  section  data. 
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Fig.  3.  1:  Flow  diagram  of  LASER  Synthesis  Code. 
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Execution  of  Analysis 

♦ 


Fig.  3.  1:  Flow  Diagram  of  LASER  Synthesis  Code  (cont'd) 
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Bogin  Execution: 
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r V tucean  ? 

Subroutine  GEAR 

Electron  cross  section  card  data  is  entered  by  DATA  BLOCK  2 to  up- 
date a main  electron  cross  section  library,  which  is  provided  on  file 
TAPE8.  The  updated  library,  which  is  generated  on  TAPE9,  may  be 
catalogued  for  future  use  if  desired.  DATA  BLOCK  2 may  be  com- 
pletely empty,  or  it  may  contain  extensive  data;  indeed,  the  initial  li- 
brary itself  can  be  constructed  by  "updating"  an  empty  file  on  TAPE8 
and  cataloging  the  resulting  TAPE9.  During  subroutine  synthesis  and 
reaction  processing,  a search  is  made  (for  each  secondary  electron 
collision  encountered)  to  determine  whether  the  appropriate  cross  sec- 
tion data  is  available  in  the  external  library.  If  it  is  not,  a warning 
diagnostic  is  generated.  The  only  cross  section  data  available  at  the 
time  of  program  synthesis  is  the  original  data  on  TAPE8.  An  omis- 
sion error,  which  would  most  likely  be  declared  fatal  initially,  could 
always  be  removed  by  an  electron  file  update  on  a subsequent  run. 

DATA  BLOCK  3 contains  comment  cards  which  are  reproduced  in  the 
output.  DATA  BLOCK  4 contains  the  experimental  input  parameters 
for  execution  of  the  analysis.  The  input  data  is  processed  and  all  of 
the  run  parameters  are  initialized  (explicitly  or  by  default).  The  mo- 
lecular kinetic  rate  constants  are  updated,  as  permitted. 

At  the  time  of  subroutine  synthesis,  reactions  are  defined  by  their 
forward  and  reverse  rate  constants.  If  a rate  constant  (e.g.,  for  the 
reverse  reaction)  was  zero,  no  translation  of  that  term  occurs  and 
no  subsequent  modification  of  that  rate  is  accessible  to  input.  For  se- 
condary electron  collisions,  however,  the  appearance  of  a zero  rate 
constant  automatically  produces  a default  to  the  Boltzmann  analysis. 

If  a nonzero  rate  constant  was  initially  provided  (at  the  time  of  synthe- 
sis) for  a secondary  electron  process,  no  coupling  of  that  process  to 
the  Boltzmann  analysis  occurs.  No  subsequent  modification  of  a rate 
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constant  which  defaults  to  the  Boltzmann  analysis  is  permitted.  An 
exception  occurs  for  the  situation  of  no  electric  discharge,  in  which 
case,  input  modification  of  a secondary  electron  rate  is  accepted.  A 
detail  balance  relation  is  enforced  for  binary  collisions,  with  the  re- 
verse rate  defined  in  terms  of  the  product  of  the  forward  rate  and  an 
exponential  factor.  Therefore,  reverse  rate  constants  for  such  pro- 
cesses cannot  be  modified  either.  An  illegal  attempt  to  modify  rates 
inaccessible  to  input  in  DATA  BLOCK  4 is  detected  (and  ignored)  and 
an  error  diagnostic  is  generated. 

The  analysis  is  executed  by  subdividing  the  total  pulse  time  TPULSE 
into  NCYCLE  (equal)  subintervals  of  time  At  = TPULSE/NC YCLE,  de- 
fining a partition  of  discrete  times  tk  = kAt  (k  = 0,  1,  2,  . . . , LIMIT). 
At  each  of  these  discrete  times,  the  Boltzmann  equation  is  solved  to 
update  the  electrical  parameters  and  electron  excitation  rates,  and  a 
variety  of  output  option  requests  can  be  specified.  Normally,  TPULSE 
is  the  total  pulse  interval,  and  LIMIT  = NCYCLE;  however,  the  com- 
bination of  parameters  TPULSE,  NCYCLE,  and  LIMIT  provides  more 
flexibility  for  the  joint  specification  of  total  integration  time  and  the 

I frequency  of  electron  kinetic  updates  and  output  generation.  Over  the 

subinterval  t^  ^ t < t^j*  the  coupled  molecular  and  circuit  equa- 
tions are  integrated  using  a multistep  Gear  technique^,  shown  in  the 
dashed  box  in  Fig.  3.1.  If  the  current  integration  step  size  is  h and 
the  time  is  t,  the  program  tests  whether  an  additional  integration  to 
time  (t  + h)  would  exceed  the  next  cycle  time  tk+1;  if  so,  program  re- 
sets the  cycle  count  to  (k  + 1)  and  extrapolates  all  parameters  from 
time  t to  time  t^^  by  the  use  of  an  array  (provided  by  Subroutine 
GEAR)  which  contains  population  densities  and  their  derivatives  (to 
various  orders).  If  (t  + h)  would  not  have  exceeded  the  next  discrete 
cycle  point,  the  GEAR  subroutine  continues  the  integration  across  the 
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subinterval.  Subroutine  GEAR  automatically  adjusts  the  step  size  h 
to  maintain  accuracy  and  stability. 

If  the  amount  of  time  consumed  by  the  Boltzmann  updates  and  output 
generation  at  a cycle  time  t^  were  negligible,  the  total  CP  time  for 
execution  of  the  complete  analysis  should  be  independent  of  NCYCLE, 
since  the  integration  would  only  be  interrupted  by  brief  pauses  at  k = 

0,  1,  2,  . . . In  actual  fact,  these  additional  tasks  performed  at  the  cy- 
cle times  t^  do  consume  only  a small  fraction  of  the  total  computa- 
tion time  for  integration  from  t^  to  tk+1.  However,  it  has  been  ob- 
served (from  actual  runs  of  the  program)  that  the  total  CP  execution 
time  seems  to  increase  proportionately  with  (rather  than  be  relative- 
ly independent  of)  NCYCLE.  At  the  discrete  time  t^,  after  the  elec- 
tron kinetic  calculations  and  output  generation,  the  GEAR  integration 
must  be  restarted  (with  JSTART  = 0).  For  the  case  of  an  electric 
discharge,  some  of  the  kinetic  rate  constants  (i.  e.  , those  for  the  sec- 
ondary electron  collisions)  change  from  t = t^  to  t = t^-,  and  it  is 
possible  that  the  restart  process  consumes  an  excessive  time  before 
the  step  size  can  be  finally  increased  to  its  most  economical  size. 

So  long  as  the  integration  loop  in  the  dashed  box  is  not  disturbed,  exe- 
cution is  efficient  and  an  economically  large  step  size  is  attained  and 
maintained.  However,  an  exit  and  subsequent  return  to  the  integra- 
tion loop  with  a restart  seems  to  cause  the  problem.  (The  problem 
does  not  seem  to  arise  for  calculations  involving  no  electric  discharge.  ) 
There  has  not  been  adequate  opportunity  to  explore  conclusively  the 
origins  of  this  difficulty,  or  to  optimize  the  execution  time.  However, 
it  is  quite  likely  that  significant  improvement  in  execution  efficiency 
should  be  possible  with  only  minimal  modifications  in  the  integration 
technique.  In  any  case,  for  the  present  version  of  the  program,  it  is 
advantageous  to  use  a relatively  small  (e.g.,  50,  100)  value  of  NCYCLE. 
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After  k = LIMIT  is  attained,  execution  terminates  and  control  is 
transfered  to  final  output  generation.  This  includes  miscellaneous 
tables  and  plots  of  electrical  and  optical  quantities,  molecular  popu- 
lation densities,  etc.  as  a function  of  time. 

The  listing  of  the  program,  which  appears  in  Vol.  II,  contains  exten- 
sive COMMENT  card  documentation  of  additional  information  which 
is  too  detailed  to  be  described  here.  As  most  of  the  FORTRAN  vari- 
able names  have  been  chosen  to  be  reasonably  descriptive  of  the  phy- 
sical quantities  they  represent,  it  is  believed  that  the  logical  and  com- 
putational structure  of  the  program  should  be  relatively  straightfor- 
ward when  supplemented  with  Fig.  3.1  and  the  present  discussion. 

The  program  has  been  written  in  such  a way  that  dimension  storage 
can  be  changed  relatively  easily.  Arrays  are  communicated  to  sub- 
routines by  variable  dimension  calls  and  thus,  it  is  sufficient  to  de- 
fine a set  of  DIMENSION  declarators  which  can  be  easily  modified. 

The  DIMENSION  statement  changes  required  are  described  by  COM- 
MENT cards  in  the  listing.  The  only  other  precaution  that  must  be 
observed  is  that  labelled  COMMON  blocks  shared  with  subroutines 
synthesized  in  previous  runs  must  be  consistent. 

External  file  usage  (TAPE2,  TAFE3,  ....  TAPE10)  is  as  follows: 
TAPE2  is  a scratch  file  used  in  subroutine  synthesis  and  analysis  of 
reaction  sensitivity.  TAPE3  and  TAPE7  are  used  for  writing  the 
FORTRAN  translation  of  the  symbolic  reaction  scheme;  they  are  sub- 
sequently combined  on  TAPE3,  which  provides  a source  file  of  80- 
BCD  character  records  of  the  statements  of  the  synthesized  subrou- 


tines JACOB,  DNDT,  and  LEVELS  which  can  be  compiled.  TAPE4 
is  a file  provided  for  writing  miscellaneous  reaction  data  associated 


with  a specified  kinetic  scheme.  TAPE4  contains  data  which  is  an 
essential  supplement  to  the  synthetically  generated  program,  and  it 
must  be  catalogued  and  saved  if  future  execution  is  intended.  TAPE3 
(which  consists  of  80- BCD  character  records  which  are  the  FORTRAN 
statements  of  the  synthesized  subroutines)  can  be  used  as  the  input 
source  for  creation  of  an  UPDATE  program  library,  which  can  be  cat- 
alogued and  saved  if  any  future  specialized  modifications  (beyond  the 
scope  of  the  synthesis  capability  of  LASER)  are  anticipated.  Normal- 
ly, it  would  be  sufficient  to  catalog  and  save  the  binary  compilation  of 
subroutines  DNDT,  JACOB,  and  LEVELS  for  future  execution.  Thus, 
even  for  specialized  situations  for  which  LASER  may  not  be  totally 
adequate  for  synthesizing  the  exact  structure  desired  for  the  rate  sub- 
routines, it  is  always  possible  to  revise  the  source  program  genera- 
ted by  means  of  subsequent  update  procedures. 

TAPE8  is  used  for  at  attached  file  of  electron  cross  section  data. 

The  electron  cross  section  library  can  be  updated  by  DATA  BLOCK  2 
card  input  to  form  a revised  file  on  TAPE9,  which  can  be  catalogued 
and  saved  for  future  use  as  desired.  TAPE10  is  the  card  input  file; 
a processing  subroutine  EDITOR  reads  card  input  data  to  an  EOF  on 
unit  10  and  generates  an  output  playback  of  the  card  images.  It  also 
copies  the  card  images  onto  TAPE5,  creating  an  effective  card  input 
file  on  unit  5.  Each  time  a new  data  block  (terminated  by  an  EOF 
card)  is  encountered,  playback  of  card  images  and  creation  of  an  ef- 
fective input  file  TAPES  occurs.  TAPE6  is  used  for  OUTPUT. 

In  the  following  sections,  some  of  the  important  subroutines  will  be 
discussed  in  more  detail.  The  card  format  for  the  input  data  decks 
will  be  described,  and  the  instructions  for  usage  illustrated  with  sam- 
ple control  cards. 
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3.2  Input  Data  Description 

Input  to  the  code  consists  of  card  and/or  tape  (disc)  data  files.  As 
the  flow  diagram  of  Fig.  3.1  indicates,  there  are  four  main  blocks 
(DATA  BLOCK  1,  2,  3,  and  4)  of  card  input;  in  addition,  an  external 
file  (TAPE8)  containing  an  electron  cross  section  library  may  also 
be  attached.  In  this  section,  a description  of  the  card  input  sequen- 
ces, with  format  specifications  and  physical  units  assumed,  will  be 
given.  The  external  cross  section  file  (attached  on  TAPE8)  is  origi- 
nally created  with  card  input  from  DATA  BLOCK  2. 


In  order  to  provide  maximum  flexibility  for  data  input  with  a minimum 
number  of  format  restrictions,  NAMELIST  entry  is  frequently  used. 
Various  constants  or  experimental  input  parameters  are  initially  as- 
signed default  values  (or  subsequently  induce  the  program  to  take  a 
suitable  default  action)  that  would  be  most  appropriate  for  the  user 
in  typical  problems  of  interest,  were  he  to  fail  to  define  certain  para- 
meters explicitly  by  input.  This  allows  for  considerable  abbreviation 
of  input,  since  only  essential  variables  need  to  be  entered  to  specify 
execution  conditions.  At  the  same  time,  however,  a wide  range  of 
variables  remain  accessible  to  input  specification  in  case  their  de- 
fault values  are  not  acceptable.  Input  format  for  NAMELIST  entry  is 
very  free,  and  eliminates  the  possibility  of  entering  numerical  data 
in  an  incorrect  field  of  columns.  Each  of  the  NAMELIST  sequences 
have  been  given  physically  descriptive  names,  and  contain  variables 
that  specify  related  parameters  for  a particular  aspect  of  the  analy- 
sis. 

In  order  to  read  data  using  NAMELIST,  input  data  must  be  provided 
in  a special  form.  The  first  character  (column  1)  on  each  card  to  be 


89 


p T 

f f 

I 


f 

\ 

I 

I 


I 


read  must  be  a blank.  The  second  character  on  the  first  card  (of  a 
group  of  data  cards)  must  be  $ (dollar  sign)  immediately  followed  by 
the  NAMELIST  name,  which  must  not  contain  any  embedded  blanks. 
The  name  must  be  followed  by  at  least  one  blank,  which  is  followed 
by  a string  of  data  items  separated  by  commas,  and  the  end  of  a data 
sequence  occurs  when  another  $ is  encountered.  Any  (or  all)  of  the 
variables  included  in  the  NAMELIST  group  can  be  entered  (without 
regard  to  order)  as  specified  input.  The  form  of  the  data  items  on  an 
input  card  is  "A  = const",  where  A is  the  FORTRAN  name  of  a varia- 
ble or  subscripted  array,  and  "const"  is  an  integer,  real,  complex, 
or  logical  conslant  (or  string  of  constants  separated  by  commas  if  A 
is  the  name  of  an  array).  Numerous  examples  will  occur  in  the  sam- 
ple output  to  be  presented  later. 

Card  input  consists  of  four  basic  blocks  of  data  to  define  the  molecu- 
lar kinetic  reaction  scheme  and  its  rate  constants,  modifications  (if 
any)  for  updating  an  external  electron  cross  section  file,  comment 
cards  to  describe  miscellaneous  information  for  reproduction  in  the 
output,  and  experimental  input  parameter  s for  execution  of  the  anal- 
ysis. 

DATA  BLOCK  J. 

DATA  BLOCK  1,  whose  structure  is  shown  in  Table  3.1,  occurs  only 
for  the  initial  synthesis  of  subroutines  and  generation  of  data  file  (on 
TAPE4)  for  a new  molecular  kinetic  reaction  scheme.  The  first  card 
is  a "header"  containing  three  (10- BCD  character)  computer  words 
which  "name"  the  new  program.  They  are  used  as  a three- line  title 
on  a cover  page  preceding  the  output  whenever  the  code  is  subsequent- 
ly executed.  Following  the  header  card,  there  is  an  arbitrarily  long  se- 
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Table  3.1:  DATA  BLOCK  1 
(Subroutine  Synthesis) 


Card  1 TITLE  (3)  (3A10) 

Header  card  contains  three  (10  BCD  char- 
acter) computer  words  which  "name"  the 
program  to  be  synthesized  from  the  reac- 
tion scheme  which  follows. 


Card  2 
Card  3 


Specification  of  reaction  No.  1. 


Card  4 
Card  5 


Specification  of  reaction  No.  2. 


Card  2n  KINETIC  (6)  (6A10) 

Symbolic  expression,  reaction  No.  n,  sub- 
ject to  certain  (minimal)  restrictions  on 
syntax,  and  miscellaneous  conventions. 

Card  2n+l  KF,  KR,  COMMENT  (5)  (2E10.3,  5X,  5A10  ) 

Forward  and  reverse  rate  constants  for  the 
reaction  (with  appropriate  units),  and  op- 
tional comment  or  reference  for  the  data. 


EOF  card  terminates  data  package.  Number  of  card 
pairs  is  arbitrary,  although  an  error  condition  results 
if  any  dimension  storage  limit  is  exceeded. 
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quence  of  pairs  of  cards  which  define  the  input  reaction  scheme,  ulti- 
mately terminated  by  an  EOF  card.  (In  practice,  the  length  of  the  re- 
action queue  is  limited  by  the  dimension  declarator  KMAX,  and  the 
number  of  new  species  encountered  during  processing  is  limited  to 
NMAX.  KMAX,  NMAX,  and  all  other  dimension  declarators  can  be 
easily  changed  in  the  main  program  LASER.  If  the  number  of  reac- 
tions or  species  exceed  the  limits,  they  are  ignored,  and  a diagnostic 
is  issued.) 

The  first  card  of  a reaction  pair  consists  of  a Hollerith  string  of  up 
to  60  BCD  characters,  containing  an  expression  of  the  form 

A + B + C + ....  r*  X + Y+  Z+  ... 

The  syntax  for  the  reaction  is  quite  flexible,  allowing  it  to  be  written 
just  as  it  would  be  in  usual  technical  notation.  Each  side  of  the  reac- 
tion may  contain  up  to  five  species  whose  names  are  defined  by  the 
first  NSIZE  (=  10)  nonblank  characters  recognized  before  a delimiter 
(+  or  r»)  is  encountered.  Embedded  blanks  are  ignored,  and  certain 
conventions  must  be  observed: 

i)  Secondary  electrons  are  denoted  by  "E",  "E-",  or  "E(-)M. 

ii)  High  energy  (e'-beam)  electrons  are  denoted  by  "HE-". 

iii)  All  other  charged  particles  must  explicitly  exhibit  their 
charge  in  their  name  (e.g.,  "F-",  "KR2(+)",  UAR(++)", 
etc.)  and  multiple  charge  must  be  explicit  (i.e.,  "AR(++)" 
rather  than  "AR(2+)"  ).  If  an  ion  is  positive,  the  "+"  sign 
must  be  immediately  followed  by  another  ,,+"  or  by  ")"  in 
order  to  avoid  confusion  with  the  normal  usage  of  "+"  in 
writing  reactions. 

iv)  Numerical  coefficient  in  front  of  the  species  not  permitted 
(e.g.  A+  B r*  C + C,  but  not  A + B r*  2C,  is  allowed). 

v)  Buffer  gases  (whose  identity  is  not  important)  can  be  de- 
noted by  "M"  (resulting  in  a total  pressure  dependence). 
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v)  Spontaneous  emission  processes  do  not  have  to  explicitly 
include  a photon  on  the  right  hand  side  of  the  reaction  un- 
less the  process  is  considered  to  contribute  "noise"  to  the 
buildup  of  an  optical  field  by  stimulated  amplification,  in 
which  case  the  photon  must  be  included  and  denoted  "HNU". 

vi)  Stimulated  emission  or  absorption  processes  are  denoted 
by  "RAD".  Except  for  the  laser  transition  involving  stim- 
ulated emission  (in  which  case  "RAD"  must  appear  on  both 
sides  of  the  reaction),  "RAD"  can  occur  only  on  the  left 
hand  side. 

As  they  are  processed,  the  reactions  of  the  input  kinetic  scheme  are 
subjected  to  several  tests  to  detect  errors.  Among  the  possible  er- 
ror conditions  encountered  are  the  following: 

i)  High  energy  electrons  must  balance  on  left  and  right. 

ii)  Buffer  gas  "M"  must  balance  on  both  sides  of  equation. 

iii)  Charge  conservation  must  not  be  violated. 

iv)  Reverse  processes  for  spontaneous  emission  not  allowed 
(i.e.,  "HNU"  must  not  appear  on  the  left  side  of  a reaction) 

v)  Stimulated  radiation  "RAD"  must  not  appear  (alone)  on 
the  right  hand  side  of  a reaction. 

vi)  Duplicate  reactions  (even  written  backward)  are  ignored. 

vii)  Detail  balance  for  binary  molecular  collisions  enforced. 

viii)  For  secondary  electron  collision  processes  which  default 
to  Boltzmann  calculation  for  rate  constant,  the  attached 
electron  cross  section  data  file  is  searched  to  verify  that 
appropriate  cross  section  data  is  available. 

ix)  Dimension  storage  (defined  by  the  declarators  KMAX, 

NMAX,  and  NKMAX)  must  not  be  exceeded  during  program 
synthesis. 

x)  Both  forward  and  reverse  rate  constants  zero;  reaction  is 
ignored  unless  it  is  a secondary  electron  collision  process, 
in  which  case  rate  constants  are  automatically  linked  to  a 
coupled  Boltzmann  analysis. 

xi)  Reaction  syntax  unrecognizable. 

xii)  More  than  five  species  on  left  or  right  hand  side. 
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If  an  error  is  encountered,  a diagnostic  warning  is  generated  in  the 

! edited  summary  of  reaction  processing,  and  the  faulty  reaction  is  ig- 

nored for  the  construction  of  the  synthetic  subroutines  and  associated 

I 

data  file.  Thus,  processing  continues  uninterrupted  throughout  the 
entire  input  reaction  queue,  and  the  final  program  generated  can  even 
be  executed  unless  the  user  declares  the  occurrence  of  such  errors 
to  be  treated  as  fatal. 

The  second  card  in  the  pair  of  cards  defining  a reaction  contains  (an 
initial  estimate  of)  the  forward  (KF)  and  reverse  (KR)  rate  constants 
for  the  process,  and  an  (optional)  50- character  field  for  comments  or 
a data  reference,  which  will  be  reproduced  in  the  output.  With  cer- 
tain exceptions  to  be  described,  both  a forward  and  reverse  rate  con- 
stant  are  typically  provided  for  each  reaction.  The  units  for  these 
rate  constants  are  determined  by  the  structure  of  the  reaction,  dnd 

e -beam  collisions;  stimulated  emission,  absorption 
Spontaneous  radiative  emission 
Binary  collisions 
Three-body  interactions 

Interaction  between  (n-  1)  species. 


are: 

2 

cm 
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cm  /a 
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cm  fa 
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cm  /s 


The  program  determines  the  form  of  the  rate  expression,  Eq.  (2), 
automatically  from  the  syntax  of  the  reaction,  subject  to  the  above 
conventions,  and  it  is  the  obligation  of  the  user  to  insure  that  input 
values  have  the  proper  units.  If  KF  = 0 (or  KR  = 0),  no  transla- 
tion of  the  forward  (or  reverse)  reaction  term  occurs  in  the  synthe- 
sized rate  subroutines  and  thus,  null  operations  (containing  unneces 


aary  multiplications  by  zero)  are  eliminated.  Rate  constants  used 
at  the  time  of  synthesis  can  be  changed  in  subsequent  execution,  how- 
ever, so  a nonzero  value  should  be  used  during  synthesis  for  any  re- 
action process  (with  the  exception  of  secondary  electron  collisions) 
which  is  not  to  be  permanently  neglected. 

For  secondary  electron  processes,  the  omission  of  input  rate  con- 
stants at  the  time  of  synthesis  automatically  results  in  the  linkage  of 
the  rate  constants  to  a coupled  Boltzmann  electron  kinetics  analysis. 
Likewise,  if  a rate  constant  is  explicitly  specified,  no  such  linkage 
occurs.  During  program  generation,  for  those  electron  collisions 
which  are  to  be  coupled  to  a Boltzmann  analysis,  a search  is  made 
to  determine  if  the  required  electron  cross  section  data  is  available 
on  the  library  file  attached  to  TAPE8.  If  the  cross  section  data  is 
available,  50- character  input  comment  is  overridden  by  a data  refer- 
ence extracted  from  the  library  file. 

For  a binary  molecular  collision  process,  the  program  automatically 
enforces  detail  balance  between  the  forward  and  reverse  rate  con- 
stants (although  with  no  provision  for  degeneracy  factors).  The  re- 
verse rate  is  defined  to  be  KR  = KF  exp  ( - AE/kT)  , where  AE  was 
given  by  Eq.  (6). 

Subsequent  modifications  of  rate  constants  at  the  time  of  execution 
are  permitted,  but  some  rate  constants  are  necessarily  inaccessible 
to  such  input.  If  no  translation  of  a reverse  process  occurred,  for 
example,  no  future  input  of  that  rate  constant  can  contribute.  Like- 
wise, if  rate  constants  automatically  default  to  Boltzmann  calcula- 
tions or  a detail  balance  expression,  opportunities  for  input  modifi- 
cation cannot  be  provided.  There  is  one  logical  exception:  if  the 
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analysis  has  been  synthesized  with  default  to  Boltzmann  electron  ki- 
netics, rate  constants  for  secondary  electron  collisions  can  be  modi- 
fied (i.e.,  defined  by  input)  whenever  the  program  is  executed  under 
conditions  of  no  electric  discharge.  If  an  illegal  attempt  is  made  to 
modify  rate  constants  which  are  inaccessible  to  input,  a warning  di- 
agnostic is  issued  and  the  illegal  attempt  is  ignored  unless  declared 
to  be  fatal. 

There  are  three  subroutines  which  are  synthesized  by  translation  of 
the  symbolic  input  reaction  scheme  into  computer- coded  equations: 

SUBROUTINE  DNDT  (N,  T,  NO,  NDOT) 

This  subroutine  calculates  the  instantaneous  rates  of  change  of  all  of 
the  dependent  variables  in  the  (first  order)  system  of  coupled  differ- 
ential equations  (consisting  of  the  master  equation  (4)  and  the  circuit 
equations  (93)  or  (94)).  Miscellaneous  optical  (gain,  absorption)  and 
electrical  parameters  (e  -beam  voltage  drop  dV/dx  and  secondary  e 

creation  rate  S,  , zero- energy  e creation  rate  S ) are  also  obtained, 
b o 

Input  Parameters: 

N = Number  of  equations 

T = Independent  variable,  time  (s) 

NO(I)  = Vector  of  dependent  variables:  1=1,  NTYPE 
are  population  densities  (cm’^);  NTYPE  +1  is 
charge,  q (Coul),  and  NTYPE+2  is  current,  i (A) 

Output  Parameters: 

NDOT  = Vector  of  rates  of  change  of  dependent  variables: 
NDOT  (I)  = dNO(I)/dt. 
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In  addition  to  the  formal  calling  parameters,  there  is  also  communi- 
cation between  DNDT  and  the  main  program  LASER  through  labelled 


COMMON  blocks.  Since  the  integration  subroutine  GEAR  is  concerned 
only  with  the  formal  calling  parameters  shown,  it  need  not  be  explicit- 
ly provided  with  miscellaneous  experimental  constants  and  input  para- 
meters available  in  the  main  program.  Likewise,  there  are  certain 
quantities  which  are  useful  by-products  of  the  computations  in  DNDT 
such  as  SB  (SJ,  SO  (S  ),  DVDX  (dV/dx),  and  RATE (I)  (the  rates  of 
Eq.  (2)).  Their  inclusion  in  labelled  COMMON  is  convenient  for  the 
direct  communication  of  information  to  the  main  program  LASER. 

SUBROUTINE  JACOB  (N,  T,  NO,  PHI ) 


This  subroutine  calculates  a Jacobian  required  by  the  integration  sub- 
routine (GEAR).  It  computes  the  (instantaneous  values)  of  the  partial 
derivatives  of  the  rates  NDOT(I)  with  respect  to  the  dependent  varia- 
bles NO(J)  at  time  T.  The  calculation  of  the  Jacobian  is  based  upon 
exact  expressions  obtained  from  formal  differentiation  of  the  rate  equa- 
tions (4)  and  (93)  (or  (94)),  rather  than  by  numerical  differencing.  Its 
availability  allows  Subroutine  GEAR  to  be  utilized  with  a "method  pa- 
rameter" of  MF  = 1.  Subroutine  JACOB  and  DNDT  share  the  same 
labelled  COMMON  blocks  with  the  main  program,  although  no  compu- 
tational results  of  JACOB  are  used  by  LASER  itself.  The  input  para- 
meters for  JACOB  are  the  same  as  for  DNDT  described  above. 

Output  Parameters: 

PHI  (I,  J)  = Jacobian  array,  dimensioned  PHL(N,  . . ) in  the 
calling  program,  defined  by 

PHI  (I,  J)  = d NDOT  (I)  / d NO  ( J) 
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SUBROUTINE  LEVELS  (Nl,  N2,  NO) 


This  subroutine  establishes  the  correspondence  between  the  lower 
and  upper  levels  of  the  inelastic  electron  collision  processes,  and 
the  species  of  the  molecular  kinetics  reaction  scheme.  It  provides 
part  of  the  linkage  which  couples  the  Boltzmann  analysis  to  the  mo- 
lecular kinetics.  Immediately  prior  to  a call  to  Subroutine  BOLTZ, 
this  subroutine  must  be  called  to  determine  the  (instantaneous)  val- 
ues of  the  lower  and  upper  state  population  densities. 

Input  Parameters: 

NO  (I)  = Vector  of  (instantaneous)  population  densities 

(cm  ) for  the  Ith  species  in  the  kinetic  scheme. 

Output  Parameters: 

Nl  (I)  = Vector  of  population  densities  for  the  lower 

states  of  the  Ith  inelastic  electron  collisions. 

I' 

N2(I)  = Vector  of  population  densities  for  the  upper 

states  of  the  Ith  inelastic  electron  collisions. 

Sample  output,  obtained  from  synthesis  of  subroutines  DNDT,  JACOB, 
and  LEVELS  for  the  KrF  laser  system,  will  be  presented  later.  The 
results  of  translation  of  the  molecular  kinetic  scheme  into  computer- 
coded  FORTRAN  equations  will  be  illustrated  for  several  reactions. 

In  addition  to  reaction  translation,  there  is  also  generation  of  inter- 
nal COMMENT  card  documentation,  which  makes  the  synthesized 
subroutines  completely  readable  and  understandable.  They  can  be 
used  as  the  source  for  creation  of  a BCD  UPDATE  program  library, 
which  can  be  catalogued  and  saved,  if  desired.  If  the  reaction  scheme 
is  to  be  subsequently  executed,  the  binary  compilation  should  be  saved. 
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DATA  BLOCK  2 


The  purpose  of  DATA  BLOCK  2 is  to  provide  updates  or  modifications 
to  the  electron  cross  section  library,  attached  on  TAPE8,  immediate- 
ly prior  to  execution.  The  updated  file,  generated  on  TAPE9,  is  then 
used  for  execution  of  the  analysis  and  can  be  catalogued  and  saved  for 
future  access,  if  desired.  Since  the  external  electron  cross  section 
library  can  be  originally  constructed  by  "updating"  an  empty  TAPE8 
file  with  the  totality  of  available  data,  the  format  for  reading  TAPE8, 
TAPE9,  and  DATA  BLOCK  2 is  identical. 

DATA  BLOCK  2,  whose  structure  is  shown  in  Table  3.2,  consists  of 
an  arbitrary  number  of  packages  of  electron  cross  section  data  for 
secondary  electron  collision  processes.  The  first  card  of  a package 
contains  a 60- BCD  character  Hollerith  field  in  which  the  symbolic  re- 
action is  specified.  The  syntax  for  reactions  is  quite  flexible,  and  is 
the  same  as  that  described  for  DATA  BLOCK  1 above.  Embedded 
blanks  are  ignored,  and  the  order  of  species  is  not  important.  Thus, 
when  a search  is  made  for  cross  section  data  for  an  electron  colli- 
sion process  encountered  in  the  molecular  kinetic  scheme,  it  is  not 
required  that  the  reaction  be  written  identically  as  it  occurs  in  the 
external  electron  data  file.  Only  the  content  of  two  reactions  com- 
pared must  be  identical;  for  example,  the  two  reactions 

KR2(+)  + E(-)  r*  KR*  + KR 

and 

E + KR2  (+)  r*KR  + KR* 

would  be  recognized  to  be  equivalent.  An  additional  convention  is 
used  for  specifying  momentum  transfer  cross  section  for  electron 
collisions  with  species  X: 
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Table  3.2:  DATA  BLOCK  2 
(Cross  Section.  Update) 


Package  1 

Package  2 

• 

• 

• 

Package  n: 

Card  1 

KINETIC  (6),  UNITS,  NPTS,  MONTH 

(6A10,  F7.3, 13,  A10)  this  card  contains  a 60- BCD  char- 
acter field  for  the  symbolic  specification  of  an  electron 
collision  process.  Cross  section  data  entered  on  subse- 
quent cards  is  assumed  to  be  expressed  in  units  of 

UNITS  x 10"  cm^.  Default  is  UNITS  = 1 if  unspeci- 

fied. NPTS  defines  subsequent  data  format. 

Card  2 

COMMENT  (6)  (6A10) 

60- BCD  character  field  for  comment  or  data  reference. 
This  information  is  reproduced  with  requested  tabular 
output  of  the  cross  section  data,  and  is  automatically  ex- 
tracted as  the  rate  reference  in  a synthesis  program. 

Card  3 

EV,  vector  of  m electron  energy  values  (eV), 

Card  4 

SIGMA,  vector  of  m corresponding  cross  section  data 
values.  Format  for  Card  3,4  is  (mEn.0),  where  m = 

NPTS  and  n =*  80/NPTS  (truncate  decimal  residue). 

• 

• 

• 

Card  i 

• 

• 

Cf.  description  of  Card  3, 

Card  i + 1 

• 

Cf.  description  of  Card  4, 

• 

• 

• 

Card  2k+l 

• 

Blank  card  terminates  cross  section  data  for  the  elec- 
tron collision  process  of  this  package. 

• 

• 

• 
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E + X r*  E + X 


Although  momentum  transfer  cross  sections  must  be  included  in  the 
electron  library  file,  the  explicit  inclusion  of  elastic  or  momentum 
transfer  collisions  in  the  molecular  kinetic  reaction  scheme  defined 
in  DATA  BLOCK  1 is  not  required. 

The  second  parameter  on  the  first  card,  UNITS,  is  a multiplicative 
factor  to  adjust  the  numerical  cross  section  values,  which  are  as- 
sumed to  be  entered  in  units  of  10  ^ cm  . (If  UNITS  is  not  speci- 
fied, it  defaults  to  1.0.  ) The  third  parameter,  NPTS,  specifies  a 
variable  FORMAT  (mEn.  0),  where  m = NPTS  and  n = 80/NPTS  (n 
is  an  integer  obtained  from  truncation  of  decimal  residue).  This 
FORMAT  specifies  the  number  of  data  fields  on  subsequent  energy 
and  cross  section  cards  to  be  read  from  the  current  package.  (If 
NPTS  ^ 0 or  if  NPTS  > 10,  program  sets  NPTS  = 10  by  default.) 

The  last  parameter,  MONTH,  is  a 10- BCD  character  specification 
of  the  date  (e.g.,  12/25/78  ) of  cross  section  submission  into  the  li- 
brary file.  If  it  is  not  provided,  it  is  automatically  generated  by 
Subroutine  UPDATE  at  the  time  it  is  added  to  the  library. 

The  second  card  of  a package  contains  a 60- BCD  character  field  for 
comment  or  data  reference.  It  is  reproduced  on  requests  for  out- 
put plots  or  tables  of  cross  sections,  and  is  automatically  extracted 
as  the  reference  for  TAPE4  data  in  the  generation  of  a synthetic  pro- 
gram. 

Following  the  first  two  cards,  there  is  an  arbitrary  number  of  card 
pairs  i and  i+1,  ultimately  terminated  by  a blank  card,  containing 
sequences  of  energies  and  corresponding  cross  section  values.  Card 
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i contains  a sequence  of  energies  (eV),  which  must  be  in  monotoni- 
cally  ascending  order,  followed  by  card  (i+  1),  which  contains  the 
corresponding  cross  section  values  (in  units  of  UNITS  x 10  ^ cm^). 

The  format  specification  is  (mEn.O),  as  described  above.  For  ex- 
ample, if  NPTS  = 8,  the  format  would  be  (8E10.0),  defining  eight 
10-column  fields  (into  which  F-numbers  could  be  entered,  if  desired). 

Ultimately,  the  collection  of  cross  section  packages  must  be  termi- 
nated by  an  EOF  card. 

DATA  BLOCK  3 

DATA  BLOCK  3,  shown  in  Fig.  3.3,  consists  of  an  arbitrarily  long 
sequence  of  comment  cards  which  will  be  reproduced  in  the  output. 

It  provides  a convenient  opportunity  for  documenting  the  output  of  an 
analysis  with  descriptive  remarks.  Format  is  (8A10). 

DATA  BLOCK  4 

DATA  BLOCK  4,  whose  structure  is  shown  in  Fig.  3.4,  defines  ex- 
perimental parameters,  numerical  control  parameters,  fatal  error 
designation,  data  initialization,  and  output  requests  for  execution  of 
the  synthesized  analysis.  It  consists  of  six  NAMELIST  data  groups, 
followed  by  an  arbitrarily  long  sequence  of  cards  which  initialize  the 
energies,  molecular  weights,  and  population  densities  of  molecular 
species.  The  NAMELIST  parameters  are  defined  in  Tables  3.5-3.  10, 
which  summarize  the  assumed  units  and  default  values.  The  NAME- 
LIST  groups  have  been  constructed  to  have  physically  descriptive 
names,  and  to  contain  logically  related  variables. 

The  format  for  the  sequence  of  species  cards  is  (A10,  3E10.3,  10X,  A10). 
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Table  3.3:  DATA  BLOCK  3 
(Descriptive  Comments) 


| 


Card  1 

COMMENT  (8) 

(8A10) 

Card  2 

• 

• 

COMMENT  (8) 

• 

• 

(8A10) 

■ 

• 

Card  n 

• 

• 

COMMENT  (8) 

• 

. 

(8  A10) 

• • 

EOF  card  termination. 
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TABLE  3.4:  DATA  BLOCK  4 
(Execution  Parameters) 


1)  $CONTROL  ...  $ 

Miscellaneous  parameters  related  to  numerical  control  of 
the  calculations--convergence  accuracy,  integration  order 
and  method,  electron  energy  resolution,  interpolation  or- 
der, fatal  declaration  of  error  conditions,  frequency  of  out- 
put generation,  etc.  (Table  3.5). 

2)  $PARAM  ...  $ 

Pressure,  temperature,  and  pulse  length  (Table  3.6). 

3)  $OPTICAL  ...  $ 

Optical  cavity  specifications- -loss  per  pass  or  threshold 
coefficient,  output  coupling  reflectivity,  length,  area,  solid 
angle  (Table  3.7). 

4)  $EBEAM  ...  $ 

External  e"-beam  current  density,  energy,  temporal  and 
energy  dependent  characteristics  (Table  3.8). 

5)  $CIRCUIT  ...  $ 

Capacitance,  inductance,  and  resistance  of  external  RLC 
driving  circuit,  initial  charging  voltage,  discharge  area 
and  electrode  separation  (Table  3.9). 

6)  $RATES  ...  $ 

(Permissible)  modifications  to  kinetic  rate  constants  for 
the  reaction  scheme  (Table  3.10). 

The  above  NAMELIST  cards  are  followed  by  (an  arbitrary  number 
of)  cards  for  initialization  of  species  parameters.  These  cards 
(which  are  ultimately  terminated  by  an  EOF  card)  are  of  the  form: 

NAME,  P,  E,  MOLWT,  OPTION  (A10,  3E10.  3,  10X,  A10) 

These  cards  define  the  energy  E (eV)  and  molecular  weight  MOLWT 
(g/mole)  of  the  species  NAME.  P is  the  concentration:  if  the  total 
pressure  is  specified  by  the  $PARAM$  card,  the  P values  represent 
relative  fractions.  If  the  total  pressure  has  not  been  previously  spe- 
cified, the  values  of  P represent  partial  pressures  (Torr),  unless 
their  values  are  > 1.0E  08,  in  which  case  they  represent  actual  pop- 
ulation densities  (cm*  ).  (Cf.  description  of  PTOT  and  ATM  in 
Table  3.6.)  If  OPTION  = 4HPLOT,  a plot  of  NAME  is  generated. 
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ie  required  accuracy  EPS  was  not  attained,  a diagnostic  is 
enerated,  and  pro  gram  terminates  unless  EXPAND  = .TRUE, 
efault  is  ITMAX  = 100. 
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TABLE  3.5s  NAMELIST  / CONTROL  / Parameters 

(Continued) 


summary  of  the  electron  distribution  function 
a.  function  of  energy  u.  (Cf.  Fig.  4.  13.  ) 
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TABLE  3.  6 : NAMELIST  / PARAM  / Parameters 
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sequent  input)  are  assumed  to  be  relative  fractions,  normalized 
to  a total  pressure  PTOT.  If  PTOT  = 0.0,  it  is  assumed  that 
input  concentrations  represent  partial  pressures  (Torr),  unless 
the  values  are  > 1.0E  08,  in  which  case  they  represent  actual 
population  densities  (cm^). 


TABLE  3.7:  NAMELIST /OPTICAL  / Parameters 


NO(l)  will  remain  zero  and  there  will  be  no  contribution  from 
the  stimulated  radiation  processes  in  the  kinetic  scheme.  (Cf. 
Eq.  (12),  (13),  (122),  and  (123).  ) 
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JBEAM  Jv(max)  A/cm  Maximum  amplitude  of  e*-beam  current  density.  Default  is 

JBEAM  = 0.0.  (REAL  variable) 
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The  general  laser  kinetics  program  LASER  requires  16  additional 
subroutines  (in  addition  to  the  three  which  are  synthesized,  JACOB, 
DNDT,  and  LEVELS).  An  independent  program,  ELECT,  has  also 
been  developed  specifically  for  Boltzmann  electron  kinetics  analysis 
and  uses  some  of  the  same  subroutines.  Complete  source  deck  list- 
ings for  both  main  programs  (LASER  and  ELECT),  as  well  as  for 
all  subroutines  used,  are  included  in  Vol.  II  of  the  present  report. 
Fig.  3.11  summarizes  the  subroutines  required  for  each  program. 
In  this  section,  some  of  the  more  important  subroutines  will  be  de- 
scribed. Program  ELECT  will  be  described  in  Sec.  3.4  below. 


SUBROUTINE  SYNTH  (LTAPE,  MTAPE,  NTAPE,  NSCRTCH,  NDATA, 
NSIZE,  MAXGAS,  GAS,  KM  AX,  NKMAX,  LEVI,  LEV2,  DATE) 


This  subroutine  translates  the  symbolic  reactions  of  DATA  BLOCK  1 
into  computer- coded  equations  for  the  synthesis  of  subroutines  DNDT, 
JACOB,  and  LEVELS  described  in  Sec.  3.2.  It  also  compiles  the 
data  defining  rate  constants  and  references  and  the  structure  of  the 
molecular  kinetic  reaction  scheme  to  generate  an  associated  file  on 
unit  NSCRTCH  ( = 4).  Source  code  is  generated  as  80- BCD  charac- 
ter records  on  unit  MTAPE  ( = 3),  which  is  compiled  for  execution. 

If  future  execution  of  a laser  kinetics  analysis  based  upon  the  input 
reaction  scheme  is  anticipated,  the  data  file  on  TAPE4  and  the  source 
code  (and/or  its  binary  compilation)  should  be  catalogued  and  saved. 

Input  and  output  parameters  for  Subroutine  SYNTH  are  summarized 
in  Table  3.12.  Specification  of  the  record  content  of  file  NSCRTCH 
is  most  easily  obtained  directly  from  the  listing. 
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Main  Program 


ANALYZE 


Synthesized  Subroutines 


DNDT 


Electron  Kinetics  Subroutines 


Utility  Subroutines 


EDITOR 


The  main  purpose  of  Subroutine  SYNTH  is  to  dissect  the  syntax  of 
(certain  types  of)  symbolic  reactions,  and  to  translate  the  fundamen- 
tal information  content  into  computer- coded  equations.  It  has  been 
"taught"  to  recognize  the  meaning  of  a reasonably  general  class  of 
expressions  as  well  as  to  detect  the  occurrence  of  miscellaneous  er- 
ror conditions.  If  necessary,  the  user  could  extend  its  translation 
or  error  detection  capabilities,  and  is  at  liberty  to  specify  in  detail 
the  structure  of  the  synthesized  subroutines.  For  example,  the  pres- 
ent techniques  could  be  extended  in  such  a way  that  an  encountered 
"species"  such  as  "XYZ(I)"  would  be  recognized  as  a representative 
member  of  a manifold  of  vibrational  levels,  and  the  reaction  syntax, 

XYZ(I)  + AB ( J)  r*  XYZ(I-l)  + AB(J+1) 

could  be  interpreted  to  be  a class  of  vibrational-vibrational  exchange 
collisions  to  be  translated  into  appropriate  DO- LOOPS  over  I,  J in 
the  generated  FORTRAN  source  code.  The  translation  of  reactions 
performed  by  the  present  version  of  SYNTH  is  best  illustrated  by  ex- 
amples and  is  therefore  postponed  until  Sec.  4,  where  results  for  the 
KrF  laser  system  are  described. 


SUBROUTINE  GEAR  (N,  T,  Y,  SAVE,  H,  HMIN,  HMAX,  EPS,  MF, 
YMAX,  ERROR,  KFLAG,  JSTART,  MAXDER,  M,  PW) 


This  subroutine  was  modified  from  "Subroutine  DIFSUB",  which  ap- 
pears in  the  work  of  Gear^  on  numerical  initial  value  problems  in 
ordinary  differential  equations.  (Cf.  Ref.  6,  pp.  155-166.)  It  is  a 
subroutine  for  integration  of  a "stiff'  system  of  coupled,  first  order, 
differential  equations  with  the  capability  for  automatic  control  of  the 
step  size  and  order  of  predictor  and  corrector.  The  replacement  of 
"MATINV"  discussed  in  the  "Note  to  Program"  (cf.  p.  158,  Ref.  6) 
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has  been  implemented  by  the  use  of  Subroutine  SIMEQ  to  avoid  the 
initial  calculation  of  a matrix  inverse,  followed  by  successive  pre- 
multiplications of  a sequence  of  vectors  with  that  inverse.  Input  and 
output  parameters  are  completely  described  by  COMMENT  cards  in 
the  listing  in  Vol.  II. 


SUBROUTINE  UPDATE  (INFILE,  NT  APE,  NSCRTCH,  LIST,  DATE) 


This  subroutine  searches  two  input  data  sources,  a file  unit  INFILE 
and/or  card  input,  to  generate  an  updated  electron  cross  section  da- 
ta library  on  unit  NTAPE.  (LASER  uses  INFILE  = 8,  NTAPE  = 9. ) 
The  updated  file  on  NTAPE  contains  all  of  the  data  of  INFILE,  modi- 
fied with  the  additions  or  revisions  defined  by  card  input.  A com- 
plete cross  section  data  package  (cf.  Table  3.2  and  the  discussion  of 
DATA  BLOCK  2)  must  be  provided  for  the  initial  addition,  or  subse- 
quent  modification,  of  any  collision  process.  NSCRTCH  is  a scratch 
file.  If  LIST  = .TRUE.,  the  contents  of  the  updated  file  NTAPE  are 
generated  in  output.  DATE  is  an  (INTEGER)  variable  containing  a 
10-BCD  Hollerith  date  (e.g.  10H  10/23/75  ). 


SUBROUTINE  PLASMA  (NDATA,  MAX,  MESH,  LHS,  RHS,  PROCESS, 
EV,  F,  G,  Q,  U0,  UM,  NTYPE,  NAME,  MISSING,  ERROR,  OUT- 
SIDE, IDEG,  OUT) 


This  subroutine  scans  the  electron  cross  section  library,  on  file  unit 
NDATA,  to  extract  the  cross  section  data  for  input  reaction  defined 
by  LHS  and  RHS,  integer  variables  generated  by  Subroutine  DEKODE 
which  uniquely  specify  the  left  and  right  hand  collision  species.  The 
reaction  is  also  described  by  PROCESS  (a  4-vector),  containing  a 40 
BCD  character  (possibly  abbreviated)  symbolic  expression  of  the  re- 
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action  (used  for  output  purposes  only).  Input  and  output  parameters 
are  fully  described  in  COMMENT  cards  in  the  listing. 


SUBROUTINE  BOLTZ  (MAX,  MESH,  NK,  GAS,  FRACT,  MIX,  NMOL, 
TMOL,  ITMAX,  TMAX,  EPS,  KAPTION,  DATE,  OUT,  EVCM,  NE, 
PROCESS,  U,  Nl,  N2,  NEL,  S,  SBEAM,  SOURCE,  X,  XQ,  QM,  F, 
G,  A,  B,  VSIG,  POWER,  PCOLL,  DISCH,  DEPOSIT,  MU,  D,  EK, 
AMPS,  UBAR,  TE,  CONVRGE,  PBAL) 


This  subroutine  implements  the  numerical  solution  of  the  Boltzmann 
equation  for  the  electron  energy  distribution,  using  the  computational 
algorithms  described  in  detail  in  Sec.  2.3.  Input  parameters  are  de- 
scribed in  Table  3.13,  and  output  parameters  in  Table  3.14.  Extensive 
internal  COMMENT  card  documentation  is  also  included  in  the  listing, 
given  in  Vol.  II. 


SUBROUTINE  DEKODE  (NAME,  IMAGE,  LHS,  RHS,  LABEL,  GAS, 
NSIZE,  NTYPE,  LONG) 


This  subroutine  dissects  the  symbolic  reaction  equations  to  determine 
the  species  on  the  left  and  right  hand  sides.  A Hollerith  specification 
of  the  reaction  is  provided  as  an  input  vector,  IMAGE  (I),  1=1,  LONG. 
Embedded  blanks  are  permitted,  and  up  to  NSIZE  characters  of  a spe- 
cies name  are  recognized  and  the  rest  ignored  (NSIZE  ^ 10).  The 
vector  NAME  provided  by  the  calling  program  contains  the  names  of 
NTYPE  species  which  have  previously  been  encountered  and  stored; 
if  additional  new  species  are  recognized  in  the  reaction  defined  by 
IMAGE,  the  vector  NAME  is  automatically  extended  to  include  them, 
and  NTYPE  is  increased  appropriately  upon  return.  GAS(J,K)  is  an 
(INTEGER)  array  containing  the  (10  BCD  character)  Hollerith  name 
of  the  Jth  species  (J  = 1,  5)  on  the  left  (K  = 1)  and  right  (K  = 2)  side 
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TABLE  3.  13s  Input  Parameters  to  Subroutine  BOLTZ 
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TABLE  3.13:  Input  Parameters  to  Subroutine  BOLTZ 

(Continued) 


DATE  INTEGER  Variable,  containing  the  date  (e.  g. , 12/15/78). 


TABLE  3.13:  Input  Parameters  to  Subroutine  BOLTZ 

(Continued) 


Secondary  electron  density. 


TABLE  3.13:  Input  Parameters  to  Subroutine  BOLTZ 

(Continued) 
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TABLE  3.14:  Output  Parameters  for  Subroutine  BOLTZ 

(Continued) 
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Electron  mobility  (cf.  Eq.  (50)). 


TABL.E3.  14:  Output  Parameters  for  Subroutine  BOLTZ 

(Continued) 


Upon  successful  convergence,  F(I)  is  the  normalized  electron 
distribution;  G(I)  = F(I)/F(1)  = fQ(u)/fo(0);  and  B(I)  = fQ(u)/fB 
(u,Te),  where  fB(u,Te)  is  the  Boltzmann  distribution  at  the  ef- 


T 

of  the  reaction,  and  LABEL  (J,  K)  are  the  integer  labels  for  GAS(J,K) 
■which  correspond  to  the  sequential  ordering  of  the  species  in  the  vec- 
tor NAME  (I),  1=1,  NTYPE.  If  there  are  less  than  the  five  species 
permitted  on  the  left  (or  right)  hand  side,  the  default  value  for  arrays 
GAS  and  LABEL  are  1H  (blank)  and  0,  respectively.  Further  details 
concerning  reaction  syntax  were  given  in  the  discussion  of  DATA 
BLOCK  1 above. 

In  order  to  make  possible  a simple  comparison  between  two  sequences 
of  collision  partners  (which  may  be  only  rearranged).  Subroutine  DE- 
KODE  also  generates  INTEGER  variables  LHS  and  RHS  to  uniquely  de- 
scribe the  left  and  right  hand  side  of  the  reaction.  If  a sequence  of  col 
lision  partners  contains  species  defined  by  the  integer  indices  two 
sums  are  defined 


LSQ  = £*/  ■ 
J 


L and  LSQ  are  then  encoded  into  an  (INTEGER)  variable,  with  format 
(14,16),  to  define  an  "identifier"  (i.e.,  LHS  or  RHS),  which  is  invar- 
iant to  a permutation  of  the  component  species  f j and  believed  to  be 
unique,  for  all  practical  purposes.  The  identifiers  LHS  and  RHSare 
used  as  Hollerith  variables,  and  contain  numbers  with  no  physical 
significance. 

LASER  initializes  the  vector  NAME (I)  in  such  a way  that  1=1  and  2 
are  reserved  for  radiation  and  secondary  electrons,  respectively. 
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Miscellaneous  Utility  Subroutines 


SUBROUTINE  SIMEQ  (A,  M,  N,  Nl,  SING) 


The  purpose  of  this  subroutine  is  to  solve  an  N x N system  of  simul- 
taneous equations  of  the  form 

^A(I,J)X(J)  = B (I), 

J 

where  I = 1,  2,  3,  ....  N.  Input  consists  of  putting  the  N x N matrix 
of  coefficients  into  the  upper  left  hand  box  of  the  array  A,  which  is 
dimensioned  A(M,NCOL)  in  the  calling  program.  Nl  different  B- 
vectors  can  be  specified,  stored  in  successive  columns  to  the  right 
of  the  N x N matrix  (cf.  Fig.  3.2).  (M  > N,  and  N + Nl  ^ NCOL. ) 
Upon  output,  the  coefficients  in  the  upper  left  hand  N x N box  are  de- 
stroyed, and  the  solution  vectors  X replace  the  corresponding  input 
vectors  B.  That  is,  the  Nl  column  vectors  B of  length  N form  an 
N x Nl  array  B,  which  is  transformed  to  A"*B  upon  output.  In  par- 
ticular, if  N B-vectors  are  chosen  to  form  the  unit  matrix,  output 
will  consist  of  the  inverse  A SING  is  a logical  variable  which  is 
returned  .FALSE,  unless  the  coefficient  matrix  is  singular. 


N f- Nl  *■ 


to##:##:# 

B 

B 

* • • 

B 

NCOL 


Fig.  3.2:  Input  array  structure  for  SIMEQ. 
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SUBROUTINE  PLOT  (MM,  MP,  MULT,  Y,  YO,  DY,  X,  XO,  DX, 
SCALEX,  SCALEY,  SAME,  CLEAR,  CENTER,  NAME,  NP,  IP) 


This  subroutine  generates  linear,  logarithmic,  or  semilogarithmic 
plots  for  NP  vectors  (NP  < 10),  Y (I,  J),  J = 1,2,  ...,NP.  A variety 
of  options  exist  for  automatic  scaling  of  the  X-  and  Y-axes  to  conve- 
nient (integer)  tick  mark  divisions,  for  simultaneous  or  independent 
scaling  of  data,  and  for  output  labeling.  Extensive  documentation, 
definition  of  input  parameters,  and  instructions  for  usage  can  be 
found  in  the  COMMENT  cards  in  the  listing. 


SUBROUTINE  AXIS  (SCALE,  XMAX,  XMIN,  XO,  DX,  XDC  ) 


This  subroutine  is  used  for  automatic  scaling  of  axes  for  Subroutine 
PLOT.  XMAX  and  XMIN  are  input  values  of  the  maximum  and  mini- 
mum values  of  the  vector  to  be  scaled.  The  range  of  values  (XMAX 
- XMIN)  is  to  be  spanned  by  a convenient  origin  XO  and  integer  tick- 
mark  spacing  DX.  XDC  is  set  equal  to  zero  unless  the  range  10  DX 
< X0/1000,  in  which  case  XDC  is  set  equal  to  XO,  XO  is  reset  to 
zero,  and  XMAX  and  XMIN  are  reduced  by  an  amount  XDC.  Thus, 
XDC  corresponds  to  a "dc  baseline"  that  is  returned  nonzero  when- 
ever the  range  of  the  plot  is  very  small  relative  to  the  absolute  mag- 
nitude of  the  values  themselves. 


SUBROUTINE  INTERP  (IDEG,  XP,  YP,  X,  Y,  MULT,  N) 


This  subroutine  will  interpolate  a vector  Y (I),  assumed  to  be  a func- 
tion of  a vector  X(I),  to  produce  that  value  YP  which  corresponds  to 
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the  value  XP.  IDEG  = 1,  2,  3,  . . . specifies  linear,  quadratic,  . . . in- 
terpolation. The  vectors  X(I)  and  Y (I)  provided  by  the  calling  pro- 
gram are  sampled  with  a repetition  index  MULT,  with  N points  de- 
fined at  I = 1,  ( 1 + MULT),  ...  ( 1 + (N-  1)  MULT)  used.  (Normal  usage 
is  MULT  =1.)  Values  of  X(I)  must  be  in  ascending  order.  XP  need 
not  lie  within  the  range  of  the  vector  X(I),  but  if  it  falls  outside,  lin- 
ear interpolation  is  automatically  used. 


SUBROUTINE  SIMPSON  (F,  M,  H,  ANS) 


This  subroutine  integrates  (Simpson  method)  a function  defined  as  a 
vector  F over  N = 2M  subintervals  of  width  H,  with  result  ANS. 


SUBROUTINE  EDITOR  (INPUT,  LIST) 


This  subroutine  reads  80- BCD  character  records  on  unit  INPUT  to  an 
EOF  and  rewrites  them  onto  TAPE5,  creating  an  effective  card  input 
file  on  TAPES  each  time  a data  block  is  processed.  LIST  = .TRUE, 
produces  an  output  "playback"  listing  of  the  input  card  images,  thus 
providing  an  exact  output  record  of  the  input  card  decks. 


SUBROUTINE  COVER  (TITLE,  NPAGE ) 


NPAGE  identical  cover  pages,  each  containing  a three-line  block-let- 
tered title  defined  by  TITLE (3),  are  generated  for  output. 


SUBROUTINE  HEADINX  (JSiMB,  JPAGE,  MESSAGE) 


This  subroutine  converts  a 10- character  phrase  defined  by  an  input 
vector  MESSAGE (10)  to  block- letters.  Refer  to  listing  for  usage. 


139 


f 


i 


; 


3.4  Boltzmann  Analysis  Program 

An  independent  program  (ELECT)  has  been  developed,  concurrently 
with  LASER,  for  Boltzmann  analysis  of  a specified  gas  mixture  and 
(sequence  of)  E/N  values.  As  Table  3.11  shows,  these  two  programs 
share  the  same  subroutine  file  for  those  computational  tasks  which 
are  common  to  both.  Furthermore,  the  structure  of  the  external  el- 
ectron cross  section  library  is  compatible  for  use  with  either  pro- 
gram (cf.  description  of  DATA  BLOCK  2 in  Sec.  3.2).  A complete 
listing  of  ELECT,  documented  by  COMMENT  cards,  is  given  in  Vol.  II. 

Input  deck  structure  for  execution  of  the  analysis  consists  of  a set  of 
data  blocks,  each  terminated  by  an  EOF  (7/8/9  card),  as  shown  in 
Fig.  3.3.  The  first  data  block  is  identical  to  DATA  BLOCK  2 (for  an 
update  of  the  electron  cross  sections)  described  in  Sec.  3.  2.  The 
second,  DATA  BLOCK  A,  consists  of  an  (optional)  sequence  of  sym- 
bolic reactions  (60  BCD  characters)  to  specify  a set  of  secondary  el- 
ectron collisions  to  which  the  Boltzmann  kinetics  analysis  will  be  re- 
stricted. If  none  are  specified,  the  program  will  automatically  use 
the  available  electron  cross  section  library  and  retain  all  collision 
processes  which  involve  any  of  the  species  which  are  subsequently 
initialized.  DATA  BLOCK  A is  followed  by  an  arbitrary  number  of 
data  blocks  containing  physical  parameters  for  execution.  The  syn- 
tax for  DATA  BLOCK  A is  the  same  as  that  described  previously  for 
the  reactions  in  DATA  BLOCK  1 and  DATA  BLOCK  2. 

The  structure  of  DATA  BLOCK  B containing  execution  parameters  is 
shown  in  Table  3.15.  NAMELIST  entry  with  suitable  default  initializa- 
tion is  used  for  input  of  numerical  control  parameters,  physical  data, 
and  initialization  of  the  components  of  the  mixture. 
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DATA  BLOCK  B 

Numerical  Control  Values 
Input  Execution  Parameters 
Initialize:  Gas  Components 


EXECUTE  ANALYSIS 


END? 


Program  Termination 


Fig.  3.3:  Input  data  deck  structure  for  Program  ELECT. 
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TABLE  3.15:  DATA  BLOCK  B 
(Boltzmann  Execution  Parameters) 


1) 

$CONTROL  ...  $ 

Miscellaneous  parameters  related  to  numerical  con- 
trol of  the  calculations- -convergence  accuracy,  ener- 
gy resolution,  interpolation  order,  output  options,  etc. 

2) 

$PARAM  ...  $ 

Pressure,  temperature,  and  E/N  values 

3) 

$SOURC E ...  $ 

Parameters  related  to  the  magnitudes  (SQ  and  S^)  and 
energy  dependence  of  an  external  ionization  source 
function  (e.g. , e -beam). 

4) 

Arbitrarily  long  package  of  species  cards  for  specifica- 
tion of  concentrations,  energies,  and  molecular  weights 

1)  NAME,  P,  E,  MOLWT  (A10,  3E10.3) 

2)  NAME,  P,  E,  MOLWT 


i)  NAME,  P,  E,  MOLWT 


These  cards  define  the  energy  E (eV),  molecular  weight 
MOLWT  (g/mole),  and  concentration  P of  the  species 
NAME.  If  the  total  pressure  has  been  specified  by  the 
$PARAM  ...  $ entry,  the  P values  represent  partial 
fractions;  if  not,  they  represent  partial  pressures (Torr) 
unless  their  values  are  > 1.0E  08,  in  which  case  they 
represent  actual  population  densities  (cm- 3). 
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3.  5 Execution  Control  Cards 


f 


In  the  present  section,  some  illustrative  examples  of  the  typical  con- 
trol card  sequences  required  for  execution  of  Program  LASER  and 
ELECT  will  be  given.  Obviously,  job  control  language  is  character- 
ized by  considerable  variations  between  different  computer  systems, 
as  well  as  by  fairly  rapid  evolutionary  changes  for  any  given  system. 
Since  the  present  programs  were  developed  for  execution  under  Ex- 
tended Fortran  Compilation  on  the  CDC  6000  and  CYBER  Series  of 
computers,  the  control  card  syntax  appropriate  for  that  system  will 
be  illustrated. 

First  of  all,  it  is  convenient  to  create,  from  the  FORTRAN  source 
deck,  a BCD  UPDATE  program  library  for  convenience  in  modifying 
the  program  if  necessary.  Fig.  3.4a  illustrates  the  creation  of  such 
a file,  which  is  given  the  name  "LASERBCD"  and  catalogued  and 
saved  for  future  access.  Furthermore,  the  binary  compilation  file 
has  been  catalogued  and  saved  as  "LASERFTN"  for  future  execution 
of  the  code.  A similar  sequence  of  operations  for  Program  ELECT 
is  shown  in  Fig.  3.4b. 

Assuming  that  the  user  finds  the  present  version  of  the  source  code 
satisfactory,  it  can  be  used  for  synthesis  and  execution  of  analysis 
with  the  control  cards  and  deck  structure  shown  in  Fig.  3.5.  The  in- 
put reaction  scheme  specified  in  DATA  BLOCK  1 is  translated  into 
FORTRAN  subroutines  on  TAPE3,  which  is  subsequently  used  as  the 
source  for  creation  of  a BCD  UPDATE  program  library  that  is  cata- 
logued and  saved  with  the  name  "KRFBCD".  The  file  of  data  associa- 
ted with  the  reaction  scheme  (created  on  TAPE4)  is  catalogued  and 
saved  with  the  name  "KRFDATA".  Two  binary  compilation  files  are 
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C SYNTH. TXXX.IOXXX.PX. 

USER ( ABC  » X YZ ) 

PROJECT (*NNN«NAME> 

REQUEST, NE*PL.*Pf. 

REQUEST. eiN.*Pr . 

. REQUEST, BINARY. *PF. 

_ REQUEST. TaPEA.'PF. 

Control  request  .t  apes, *pe  . 

Cards  < ATTACH (LGO.LASCRFTN, IO=XXX > 

J LOSET(PRESET=INUEF) 

LGO. 

UPDATE (N . I -T  APE3 ) 

CAT ALOG (NE«PL .KRFBCD.  tD=XX*> 
FTN(I,B=BIN,R=3,EL=F) 

CATALOG (BIN.KRF8 IN. IO=X XX) 

COPYL (LGO. B IN, B INARY,, AR) 

CATALOG (B I NARY. KRFLASER.ID=*XX) 
CATALOGCTAPEA.KRFDATA. !D=X*A) 
LDSET(PRESET=INDEF) 

BINARY <PL*25000) 

V CATAL0G(TAPE9, CROSSSECTIONS. ID=XXX) 
7/8/9  CARD 
HEADER  CARD 
A » B <•  C ♦ 0 

KF  KR  REFERENCE 


DATA 
BLOCK  1 


ARBITRARILY  LONG  SEQUENCE  of  reactions 


7/8/9  CARD 


DATA 
BLOCK  2 


STACKED  E-  CROSS  SECTION  PACKAGES 
(ARBITRARILY  MANY  ARE  PERMITTED.) 


DATA 
BLOCK  3 


7/8/9  CARO 

COMMENT  CAROS 


DATA 
BLOCK  4 


7/8/9  CARO 
SCONTROL  ...  * 
SPAR AM  ...  $ 
SOPTICAL  ...  * 
SEBEAM  ...  S 


6/7/S/9  CARO 


Fig.  3.5:  Synthesis  of  program  defined  by  input  reaction  scheme, 
catalogue  of  permanent  files,  and  execution  of  analysis. 
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created.  The  compilation  of  the  synthesized  subroutines  (i.e.  JACOB, 
DNDT,  and  LEVELS)  for  the  input  reaction  scheme  is  catalogued  as 
"KRFBIN".  The  complete  binary  file  (containing  the  main  program 
LASER,  all  of  its  subroutines,  and  the  subroutines  just  synthesized) 
is  catalogued  as  "KRFLASER".  The  analysis  which  follows  (defined 
by  DATA  BLOCK  2,  3,  4)  is  then  executed  under  the  binary  file  KRF- 
LASER. Note  that  no  electron  cross  section  library  has  been  attached 
on  TAPE8;  the  first  time  the  program  is  run,  DATA  BLOCK  2 must 
contain  card  input  of  all  available  electron  cross  section  data  for  the 
initial  creation  of  such  a library  on  TAPE9,  which  is  shown  catalogued 
as  CROSSSECTIONS.  For  all  future  execution  of  an  analysis  based 
upon  the  input  reaction  scheme,  the  binary  file  KRFLASER  can  be 
used,  and  CROSSSECTIONS  can  be  attached  to  TAPE8,  as  shown  in 
Fig.  3.6. 


If  it  is  desired  to  modify  the  present  version  of  the  program,  this  can 
easily  be  done  using  UPDATE  procedures.  If  these  modifications  in- 
clude (but  are  not  necessarily  limited  to)  changes  in  subroutine  SYNTH, 
the  control  cards  and  deck  structure  shown  in  Fig.  3.7  are  appropriate. 
The  three  previously  created  permanent  files  LASERBCD,  LASERFTN, 
and  CROSSSECTIONS  are  attached  as  designated.  Modified  subrou- 
tines are  compiled  and  the  old  versions  replaced  by  the  new  using  the 
COPYL  command.  The  new  binary  compilation  file  is  catalogued  (as 
LASERFTN)  and  the  old  one  purged.  Program  LASER  is  then  execu- 
ted in  its  updated  version,  according  to  the  input  of  DATA  BLOCK  1, 

2,  3,  and  4.  Revisions  or  additions  to  the  electron  cross  section  data 
of  DATA  BLOCK  2 are  incorporated  permanently  by  cataloging  the 
updated  file  on  TAPE9  as  "CROSSSECTIONS"  and  purging  the  initial 
file  that  was  attached  to  TAPE8.  If  modifications  to  program  LASER 
do  not  involve  subroutine  SYNTH,  and  it  is  only  desired  to  make  cer- 


Control 

Cards 


DATA 
BLOCK  2 


DATA 
BLOCK  3 


DATA 
BLOCK  4 •< 


1 

{ 


LASER. TXXX.IOXXX.PX, 

USER ( ABC. XYZ) 

project  c *nnn*name * 

ATTACH (BIN. LASERETN, ID.XXX » 
ATTACHIBINART.KREBIN. io*xxx> 

ATTACH (T AREA. KRFOAT A, IO*XXX I 
ATTACH ( T APE8.CH0SSSECT IONS . IO»XX*) 

LOAO(HINARYI 

LOSE  T ( PRESET  » INOEF  > 

BIN(PL>2S000) 

7/8/9  CARO 

stackeo  e-  cross  section  packages 

(ARBITRARILY  MANY  ARE  ^ERMlTTEO.) 


7/8/9  CARO 

COMMENT  CAROS 


7/8/9  CARO 
SCONTROL  .. 
SPARAN  ... 
(OPTICAL  .. 
(EBEAM  ... 
(CIRCUIT  .. 
(RATES  ... 


7/8/9  CARO 
6Z7/8/9  CARO 


• 

Alternate 

Control  " 

Cards 

LASER. TXXX.IOXXX.PX, 

USER ( ABC .XYZ) 

PROJECT (•NNN«NAM£I 

ATTACH(BINARY. KRFLASER, TO»**X 1 
attacm(tape*.krfoata. in.xxxi 

ATT ACH( TAPES. CROSSSECT »OHS.lO«*XX» 
LOSET (PRESET»INOEn 

BINARY (PL«2S000) 

7/8/9  CARO 

• 

4 

4 

4 

6/T/0/9  CARD 

Fig.  3.  6:  Execution  of  the  analysis  by  loading  LASERFTN  and 
KRFBIN,  or  the  complete  binary  file  KRFLASER. 
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Control 

Cards 


Update 

Program 

Library 

DATA 
BLOCK  1 

DATA 
BLOCK  2 


DATA 
BLOCK  3 


DATA 
BLOCK  4 


SYNTH. TXXX.IOXAX.PX, 

USER ( ABC.XYZ) 

PROJECT (*NNN»NAME) 

ATTACH «OLOPLtLASERBCO.TO«X**» 

ATTACH (OB J.LASERFTN.ID«XXX» 

ATT ACH< TAPES. CROSSSECT  IONS* I0»XXX) 
REQUEST, TAPE9.«PF. 

REOUEST .BIN, *PF t 
RE0UEST.8IN1 »»PF. 

REOUEST.  TaPE<*»*PF. 
request. ne*pl.*pf. 

UPDATE  (P»L*1 > 

FTNI I«C0HPILE»R«3.A.EL=F) 

RENINDI08J.L00) 

COPYL (0BJ.LG0.BIN) 
CATALOG(BIN,LASERFTN,Io*XXX) 

PURGE (OBJ) 

LDSET (PRESET* INOEF) 

BIN 

UPDATE  CN. I*TAPE3) 
CATALOG(NE«PL.KRFBCD.IO*XXX> 
FTN(t*B*BINI»R*3*EL*F) 
C*TAL0GIBIN1.KRFBIN.ID«XXX> 

COPYL (BIN.BINl .BINARY , . ARI 
CATAL0GITAPE<.»KRFDATA.ID«XX*> 

LDSET (PRESET* INOEF) 

BINARY (PL*25Q00) 

CAT AL0GCTAPE9. CROSSSECT  IONS. ID«XXX) 
PURGE (TAPES) 

7/8/9  CARO 

INSERTIONS  AND  deletions 


7/8/9  CARO 

header  card 

arbitrarily  LONG  SEQUENCE  of  reactions 


7/8/9  CARD 

stacked  e-  cross  section  packages 

(ARBITRARILY  MANY  ARE  PERMITTED) 


7/8/9  CARD 

COMMENT  CAROS 


7/8/9  CARD 

EXECUTION  DATA  PACKAGE 


7/8/9  CARO 
b/7/8/9  CARD 


Fig.  3.7:  Synthesis  and  catalog  of  program  KRFBCD  under  a 
modified  version  of  LASERBCD,  followed  by  execu- 
tion of  the  analysis. 
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— 

tain  changes  related  to  execution  of  the  analysis  (not  involving  the  syn- 
thesis of  subroutines),  Fig.  3.8  contains  the  appropriate  control  card 
and  input  deck  structure. 

It  is  conceivable  that  the  user  may  desire  to  make  changes  both  in  the 
main  program  LASER  (or  in  any  of  its  subroutines)  and  in  the  synthe- 
tic subroutines  (JACOB,  DNDT,  or  LEVELS).  To  facilitate  this,  it 
is  convenient  to  have  catalogued  the  BCD  UPDATE  program  library 
KRFBCD  (as  described  in  Fig.  3.5).  Both  LASER  BCD  and  KRFBCD 
can  be  updated,  with  their  modified  subroutines  replaced  (by  COPYL) 
on  the  binary  compilation  file  KRFLASER  used  for  execution.  It  would 
be  much  more  difficult  to  extend  the  translation  capabilities  of  subrou- 
tine SYNTH  to  automatically  recognize  new  reaction  syntax,  or  to  in- 
corporate new  computational  features  into  the  molecular  kinetic  sub- 
routines, than  it  would  be  to  merely  modify  the  synthetic  subroutines 
generated  by  the  present  program.  There  are,  therefore,  situations 
where  an  update  of  the  program  library  KRFBCD  (supplemented  by  any 
changes  required  for  insuring  compatibility  with  LASER  or  its  subrou- 
tines) would  probably  be  more  efficient  than  an  attempt  to  extend  sub- 
routine SYNTH.  On  the  other  hand,  extensions  which  are  likely  to  be 
of  general  applicability  to  future  usage  of  the  code  should  probably  be 
incorporated  directly  into  subroutine  SYNTH. 

Execution  of  the  Boltzmann  analysis  program  ELECT  can  be  carried 
out  with  the  control  card  sequence  and  data  deck  structure  shown  in 
Fig.  3.9.  The  electron  cross  section  data  library  CROSSSECTIONS 
(discussed  above)  is  compatible  with  either  program  (LASER,  ELECT) 
and  the  structure  of  the  electron  update  input  (DATA  BLOCK  2)  is  iden- 
tical for  both  programs.  Obviously,  the  same  procedures  for  program 
modification  by  update  of  BOLTZBCD  (as  discussed  for  LASERBCD  ) 
are  possible. 
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Control 

Cards 


Update 

Program 

Library 

DATA 
BLOCK  2 


LASER «TXXX» IOXXX.PX. 

USER ( ABC.XYZ) 

project («nnn«name) 

AFTACH(OLOPL*LASERBCD. tO»XAXj 
ATTACH(OBJ*LASEHFTN» IOzXXX) 

ATTACH (TAPES. CROSSSECT  TONS* In*XxX) 

ATTACH(TAPE**.KRFDATA.ln=XX*> 

ATT  ACH (b INARY .KRFBIN t IOzXX A ) 
REQUEST. TAPE9.»PF. 

REQUEST. BIN. «PF. 

UPDATE (P«LS1 1 

FTN( I«COMPILE.R=3.A.EL=F) 

REDUCE . 

REWIND(OBJ.LGO) 

COPYL (OBJ.LGO.BIN) 
catalogibin.laserftn.io=xxa) 

PURGE (OBJ) 

LOAD(BINARY) 

LDSET (PRESET- INOEE) 

BIN(PL*25000I 

CATALOG ( T APE9. CROSSSECT IOnS.jo«XXX) 

PURGE (TAPES) 

7/8/9  CARO 
•IOENT. -CHANGE- 

INSERT  IONS  AND  deletions 


7/8/9  CARD 

STACKED  E-  CROSS  SFCTlON  PACKAGES 
(ARBITRARILY  MANY  PERMITTED* 


DATA 
BLOCK  3 


7/8/9  CARD 

COMMENT  CAROS 


DATA 
BLOCK  4 


7/8/9  CARO 
(CONTROL  ...  * 
SPARAM  ...  S 
SOPTICAL  ...  * 
SEBEAM  ...  S 
(CIRCUIT  ...  * 
(RATES  ...  ( 
GAS1  PI 

GAS2  P2 

GAS3  P3 


PI 

El 

HI 

P2 

E2 

R2 

P3 

E3 

M3 

• 

• 

• 

PK 

EK 

mk 

7/8/9  CARO 
fc/7/8/9  CARO 


Fig.  3.8:  Program  library  LASERBCD  is  modified  by  update, 
previous  binary  file  KRFBIN  is  loaded,  and  analy- 
sis is  executed  under  new  version  of  LASERFTN. 
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Control 

Cards 


BOLTZ.TXXX* IOXXX.PX. 

USER ( ABC.XYZ) 

PROJECT  I »NNN*N AM£ > 

ATTACH (TAPES. CHOSSSECT IONS* IO«XX*» 
ATTaCH(LGO.BOLTZFTN, IO«XXX» 

REDUCE . 

LOSET (PRESET*INOEF> 

LGO (PL*25000 ) 

7/8/9  CARD 


DATA 
BLOCK  2 


ELECTRON  CROSS  SECTION  PACKAGES 
(ARBITRARILY  MANY  ARE  PERMITTED. t 


7/8/9  CARO 
AR  * E •>  AR*  • E 
AR  * E » AR**  ♦ E 
AR  * E » AR(»>  • E » E 


DATA 
BLOCK  A 


SECONDARY  ELECTRON  COLLISION  PROCESSES 
INCLUOEU  IN  THE  BOLTZMANN  ANALYSIS 


DATA 
BLOCK  B1 


7/8/9  CARD 
SCONTROL  ...  * 
SPAR AM  ...  S 
S SOURCE  ...  S 
GAS  1 PI  El 

GASZ  P2  E2 

GAS3  P3  E3 


DATA 
BLOCK  B2 


7/8/9  CARO 
SCONTROL  ...  S 

sparam  ...  s 

SSOURCE  ...  s 
GAS l PI  El 

GAS2  P2  E2 


7/8/9  CARO 


6/7/S/9  CARO 
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4.0  KrF  LASER  KINETICS:  Illustrative  Example  of  Program 

Synthesis  and  Execution 

The  Computer  Program  LASER  that  was  described  in  Sec.  3 can  be 
used  to  (automatically)  synthesize  and  execute  (for  a reasonably  ar- 
bitrary reaction  scheme)  a coupled  analysis  of  molecular  kinetics,  el- 
ectron kinetics,  radiative  extraction,  and  external  (RLC)  driving  cir- 
cuit. In  this  section,  usage  of  the  code  will  be  illustrated  with  sam- 
ple output  obtained  from  analysis  of  the  KrF  laser  kinetic  system1"5. 
Names  of  permanent  files  in  the  discussion  of  control  cards  in  Sec. 

3.5  were  chosen  in  anticipation  of  this  example. 

4.  1 Molecular  Kinetic  Subroutine  Generation 


As  many  as  80  reactions  have  been  proposed  as  being  of  possible  im- 
portance for  the  KrF  laser  kinetic  reaction  scheme.  Fig.  4.1  shows 
the  "playback"  (generated  in  output  by  subroutine  EDITOR)  of  the  ex- 
ac^  input  card  images  of  DATA  BLOCK  1 (used  in  the  deck  structure 
of  Fig.  3.5,  for  example)  for  synthesis  of  a KrF  laser  kinetics  analy- 
s*8,  The  header  Card  1 contains  three  (10  BCD  character)  computer 
words  ("  NORTHROP  ",  "KRF  LASER  ",  and  "KINETICS  ")  which 
name  the  generated  program  and  are  used  for  a three-line  block-let- 
ter title  in  subsequent  output.  Following  the  header  card,  there  are 
80  pairs  of  cards  which  define  the  kinetic  reaction  scheme  and  its  rate 
constants  (cf.  description  of  DATA  BLOCK  1 in  Sec.  3.2). 

The  reactions  of  Fig.  4.1  are  processed  by  Subroutine  SYNTH  and  are 
translated  into  FORTRAN  subroutines  (JACOB,  DNDT,  and  LEVELS). 
During  the  translation  process.  Subroutine  SYNTH  also  generates  its 
own  COMMENT  card  documentation  of  the  synthesized  subroutines  in 
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Fig.  4.1:  DATA  BLOCK  1 card  input  for  definition  of  KrF  kinetic  reaction  scheme. 
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order  to  make  them  completely  readable  and  accessible  to  modifica- 
tion. Fig.  4.2,  4.3,  and  4.4  present  samples  of  the  FORTRAN  state- 
ments generated  for  subroutines  DNDT,  JACOB,  and  LEVELS,  re- 
spectively. Although  it  would  be  impractical  to  reproduce  the  entire 
listings  of  these  synthesized  subroutines  here,  it  is  informative  to 
present  excerpts  which  illustrate  the  translation  for  miscellaneous 
representative  types  of  collisions,  as  is  done  in  Fig.  4.2  and  4.3. 

As  the  reaction  queue  of  DATA  BLOCK  1 is  scanned  by  Subroutine 
SYNTH,  input  reference  data,  informative  comments,  and/or  warn- 
ing diagnostics  are  generated  to  produce  an  output  summary,  shown 
in  Fig.  4.5.  A reaction  with  unacceptable  syntax  or  content  is  ignored, 
although  processing  continues.  (An  igno  'ed  reaction  is  not  assigned 
a number  in  the  summary  of  Fig.  4.5.)  After  the  entire  reaction  list 
has  been  processed,  a cross-reference  table,  illustrated  by  Fig.  4.6, 
is  generated.  This  table  provides  for  rapid  identification  of  all  re- 
actions that  involve  any  particular  species,  and  is  convenient  when- 
ever it  is  desired  to  delete  some  species  from  the  kinetic  scheme. 

After  generation  of  the  FORTRAN  source  code  on  TAPE3  and  construc- 
tion of  the  data  file  on  TAPE4,  there  is  an  exit  from  program  LASER. 

TAPE3  is  used  as  the  source  for  creation  of  a BCD  UPDATE  file  (cf. 

Fig.  3.5)  which  is  catalogued  as  KRFBCD,  and  its  binary  compilation 
is  combined  with  LASERFTN  to  create  file  KRFLASER  for  execution. 

The  data  file  on  TAPE4  is  catalogued  as  KRFDATA. 

4.2  Execution  of  the  Analysis 

With  the  use  of  permanent  files  KRFDATA  and  KRFLASER,  the  analy- 
sis can  be  executed,  either  as  a continuation  from  the  initial  synthesis 
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SUBROUTINE  Or.vT  (N,  T.  NO.  NOOT) 
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THIS  SUBROUTINE  was  synthesized  by  EDITING  AN  INPUT  FILE  of  SYM- 
BOLIC reactions  which  oefinE  a coupleo  system  of  electron  ano 
MOLECULAR  KINETICS.  IT  RETURNS  the  RATES  NOOTdl  « (0/OT  ) NO  < I ) * 

I * 1 *?* • • .NTYPE  ICH-3/SEC).  RATE  CONSTANTS  KF  ANO  KR  HAVE  UNITS 
OF  CM2.  SEC-1.  CM3/SEC.  CM6/SEC*  ...  AS  APPROPRIATE. 

THE  GENERAL  KINETICS  SYNTHESIS  PROGRAM  WHICH  AUTOMATICALLY  GEN- 

e«*teo  THir  subroutine  was  oevelopeo  by  — 


I DR.  WILLIAM  B.  LACINA 

I NORTHROP  RESEARCH  ANO  TECHNOLOGY 

I ONE  RESEARCH  PArk 

I PALOS  VEROES  PENINSULA.  CA  902TA 

I TEL*  *213)  377-4811*  EXT.  322 


DIMENSION  No II).  NOOT Cl) 

REAL  no.  NTOT.  NDOT » NOISE.  NE.  KF,  KR,  KB,  KT.  MU,  LO.  IBEAM, 

1 jbeam,  length 

COMMON  / data  / RATE!  200).  KF C 2<>0).  KRI  200),  VSIG<2.  25).  £(30) 
COMMON  f CONST  / NTOT,  TMOL » FREQ.  HNU 
COMMON  / OISCH  / L0,  CO,  R0.  Mij,  AREA.  D 

COMMON  / SOURCE  / UPLUS.  JBEAM.  OVOX.  DEPOSIT,  ENERGY,  SB.  SO 
COMMON  / GAINS  / ALPHA.  GAMMA,  GAIN,  ABSORB.  OMEGA,  LFNGTH.  CAVITY 

DATA  KB»  E0,  h.  C.  PI  / B.614E-05.  1.602E-19.  6.62SE-34. 

1 3.000E.10,  3.14159  / 

IBEaM  3 JBEAM«OEPOSlT*SHAPE*T| 

KT  3 KB#TMOl 
DO  I I * 1»N 
1 NOOTU)  * 0. 

NOOT II)  * - C*GAMMa»NO(1) 
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G*lN  * SIGMA«(N2-N1)  is  the  laser  transition  gain 

ABSORB  * SUMKJ SIGMA  IK )»NK1  IS  THE  TOTAL  ABSORPTION  OF  THE  MEDIUM 

ALPHA  * IbAIN-ABSORB)  IS  THE  NET  GAIN  IN  THE  MEDIUM 

gamma  « threshhold  gain  coefficient  icm-h 

CAVITY  * MIRROR  SEPARATION  (CM) 

LENGTH  « length  of  ACTIVE  MEOIUM  (CM) 

OMEGA  * AR£A/CAVITY»»2 

area  * AREA  of  optics  (Cm2) 
gamma  * ILOSS  » LN ( I /R ) /2 1 /LENGTH 
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Fig.  4.2:  (Synthesized)  Subroutine  DNDT  listing. 
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0OLIC  REACTIONS  WHICH  DEFINE  A COUplEO  SYSTEM  OF  ELECTRON  ANO 

JACOB 

7 

c 

MOLECULAR  KINETICS  EQUATIONS.  IT  RETURNS  THE  JACOBIAN  MATRIX. 
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PHI(I.J)  » DlNDOT(I)  )/0(N0(J)  1 

JACOB 

10 

c 

JACOB 

11 

c 

WHERE  I.J  * 1.2.3... .NT yPE.  N IS  THE  DIMENSION  DECLARATOR  FOR  PHI 
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(as  indicated  in  Fig. 3.5)  or  directly  (with  the  control  card  deck  given 
in  Fig.  3.6).  Fig.  4.7  shows  the  "playback"  (generated  in  output  by 
subroutine  EDITOR)  of  the  exact  input  card  images  of  DATA  BLOCK 
4,  which  specifies  the  experimental  conditions.  The  structure  of  the 
input  deck  has  been  described  in  Sec.  3.2  and  summarized  in  Tables 
3.4-  3.10. 

At  the  beginning  of  the  output  after  execution  of  the  analysis,  a cover 
page  as  shown  in  Fig.  4.8  (actually,  two  identical  such  pages)  is  gen- 
erated, containing  a block- letter  title  defined  by  the  header  card  of 
DATA  BLOCK  1 at  the  time  of  program  synthesis  (cf.  Fig.  4.1).  As 
described  in  Table  3.5,  the  IO  input  vector  provided  by  $CONTROL$ 
entry  (cf.  Fig.  4.7)  is  used  to  specify  output  option  requests.  IflO(lO) 
is  nonzero,  an  initial  summary  of  the  reaction  scheme  and  its  rate 
constants  (as  modified  by  $RATES$  input  entry)  is  generated,  shown 
in  Fig.  4.9.  (The  effect  of  rate  modifications  will  be  illustrated  in 
Sec.  4.3  below.  ) Only  the  first  (50  BCD  character)  line  of  comments 
which  appeared  in  Fig.  4.5  is  retained.  If  10(8)  and  10(9)  had  been 
entered  nonzero,  a tabular  summary  and  plot  of  all  inelastic  and  mo- 
mentum transfer  electron  cross  section  data  (used  in  the  analysis)  is 
provided,  as  shown  in  Fig.  4.10  and  4.11.  The  parameters  10(1),  I = 

1,  2,  . . . , 7 are  used  to  specify  the  frequency  of  certain  output  options 

described  in  Table  3.5.  Fig.  4.12-  4.  18  are  samples  of  such  output 

_q 

corresponding  to  cycle  k = 20  (t  = 200  x 10  s).  The  first  five  of 

these  figures  summarize  results  of  the  electron  kinetics  analysis  at 
t = 200  ns,  while  the  sixth  (Fig.  4.17)  summarizes  the  molecular  kine- 
tic, electrical,  and  optical  parameters. 

Note  that  Fig.  4.17  contains  a tabulation  of  the  population  densities  and 
their  instantaneous  rates  of  change  for  all  of  the  species.  Also  shown 


183 


is  an  effective  (instantaneous)  time  constant  Tj  for  each  species  I, 

Tjd)"1  = | dNjW/dtl/Njft), 

which  is  useful  for  indication  of  how  rapidly  any  given  species  is  evolv- 
ing at  any  given  time.  Fig.  4.18  is  a (partial)  reproduction  of  output 
which  summarizes  the  numerical  contribution  of  each  kinetic  reaction 
to  every  species,  and  is  useful  for  assessing  reaction  sensitivities.  If 
a particular  reaction  makes  only  a negligible  contribution  under  typical 
conditions  of  interest,  it  can  easily  be  deleted  in  a future  version  of  the 
program.  Since  the  output  format  of  Fig.  4.12-  4.18  is  essentially  self- 
explanatory,  no  further  description  should  be  necessary. 

At  the  conclusion  of  the  analysis,  numerous  parameters  are  plotted  or 
tabulated  as  a function  of  time,  such  as  e -beam  current  density  (Fig. 
4.19),  intracavity  radiation  intensity  (Fig.  4.20),  circuit  voltages  and 
power  densities  (Fig.  4.21),  instantaneous  net  gain  coefficient  (Fig.  4.22) 
and  medium  gain  and  absorption  coefficients  (Fig.  4.23),  optical  and  el- 
ectrical power  densities  (Fig.  4.24)  and  energy  densities  (Fig. 4. 25),  cir- 
cuit voltages  (Fig.  4.  26),  plasma  conductivity  and  discharge  impedance 
(Fig.  4.27),  discharge  current  density  (Fig.  4.28),  optically  extracted 
power  (Fig.  4.29)  and  energy  (Fig.  4.30)  density  and  efficiency,  and  the 
population  densities  (and  rates  of  change)  for  miscellaneous  species  as 
requested  (Fig.  4.31  - 4.35). 
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SUMMARY  Or  C*RO  IMAGES  FOR  INPUT  DATA  DECK 
IOATE>  01/12/79) 


1 2 3 * 5 6 7 ft 

1 23*567090 I 23*567890 123*567890 123*567890 1 23A567890 1234567890 123*567890 123*567890 


SCONTROL 


10 4 1 > « 10*?0*20*20« 10. 10,20,0,0, 1 • 


2 

• ••• 

ERA*  * 20.» 

3 

• • • . 

NCTCLE  * 100S 

* 

....  SpARAM 

TPULSE  * 1.0E-06. 

5 

• ••• 

ATM  « 2. OS 

6 

*•*.  SoPTlCAL 

REFLECT  * 70.. 

7 

• ••• 

LOSS  ■ 0.. 

A 

• ••• 

AREA  ■ 25.» 

9 

• ••• 

CAVITY  . 130.. 

10 

• ••• 

LENGTH  m 75. OS 

11 

....  SeBEAM 

JBEAH  > 10.0, 

12 

• • • • 

energy  « 300.. 

13 

• ••• 

UB  * 20., 

1* 

• • • • 

TB<1>  * 0.«  50.,  100..  ISO., 

200*«  250. • 300*9 

15 

• ••• 

350..  *00..  *50.. 

500**  550. t 600*9  650*9 

16 

TOO.,  750..  800., 

650*9  900.*  950*9  1000*9 

17 

• ••• 

UNITS  * 1.0E-09. 

IB 

• ••• 

JB ( 1 ) * O.f  .82*  lot  let  o^t 

•659  *77*  .73«  .69* 

19 

• ••• 

.67*  #69*  #7?t  #75»  0 75*  o7?*  #65*  #62*  •?* 

20 

.05,  .02,  0., 

21 

• • • • 

FACTOR  . 2.01  , 

22 

....  SCIWCUIT 

CAPAC  . 1.00E-06, 

23 

• ••• 

INDUCT  > 200.E-09, 

2* 

• ••• 

RESIST  * 0., 

25 

• • • . 

KVOLT  * so.. 

26 

• ••• 

AREA  * 2000., 

27 

• ••• 

OIST  • 10. S 

2B 

....  SRATES  s 

29 

.... AR 

89.9  *0. 

30 

... .KR 

10.0  a*. 

31 

••*.12 

0.1  38. 

PLOT 

32 

.... AR* 

11.5 

PLOT 

33 

... .KR* 

9.9 

PLOT 

3* 

... . AR£( * ) 

13.0 

PLOT 

35 

....KR2<.> 

11.8 

PLOT 

36 

.... arkR 1*1 

12.* 

37 

. ...AR** 

13.5 

3D 

....aR(*> 

15.8 

PLOT 

39 

... .KR** 

11.5 

*0 

,...KR|*> 

1*.* 

PLOT 

123*567890 123*567890 1 23*567890 123*567890 123*567890 1 23*567890 123*567890 I23*567»"0 

1 

2 3 * 

5 6 7 8 

SUMMARY  of  CARO  IMAGES  FOR  INPUT  OATa  OECK 

10ATEI  01/12/79) 

CARO  1 

2 3 * 

5 6 7 8 

NO 

123*567090 1 23*567890123*567890 1 23*567890 1 23*567890 1 23*567890 1 23*567890 1 23*567890 

*1 

... .ARF* 

6.5 

PLOT 

*2 

...  .krf* 

5.0 

PLOT 

*3 

• • • . AR2* 

9.5 

** 

... .KR?* 

8.2 

*5 

...,ar ar* 

5.0 

*6 

....arkRf* 

*.0 

*7 

... .KRzE* 

3.0 

*D 

....ARkR* 

8.8 

*9 

.,..f- 

PLOT 

50 

....EI-* 

PLOT 

123*567890123*567890123*567890123*567890123*567890123*567890123*567890123*567890 

I 

2 3 * 

5 6 7 8 

Fig.  4.7:  DATA  BLOCK  4 card  input  of  numerical  control 
values,  experimental  parameters,  modification 
of  rate  constants,  and  initialization  of  species. 
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Fig.  4.  8:  Cover  page,  generated  at  beginning  of  each  output,  with  a block- letter  title 
defined  by  the  header  card  of  DATA  BLOCK  1 at  time  of  program  synthesis. 
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Fig.  4.9:  Continued 
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Fig.  4.10:  Tabular  summary  of  electron 
cross  section  data. 
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Fig.  4.11:  Plot  of  electron  cross  section  data. 
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LOG  Of  ELECTRON  ENERGY  U IEVI 


TIHC  T ■ MflCI'  SEC 
GAS  MIX  TUBE  — 

All  /KM  /ri  ■ 89.90  / 18.80  / .18 

PLASMA  PARAMETERS 


E iriELOl 

■ 

8.3I23E  *02 

VOLT/CM 

TMOL 

• 

308. 

0E6  K 

NMOL 

• 

6.8B93EM9 

CM- 3 

PTOT 

• 

1S20.0 

torr 

• 

2.80 

ATM 

E/N 

■ 

l.T00ie-IT 

VOLT  CM2 

E/P 

■ 

.55 

V/CM/IORR 

• 

.A2 

kv/cm/atm 

<U> 
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3.087 

EV 

TE  • 2<U»/3M 

• 

23889. 

OCG  R 

oir fusion  o 

• 

2I31.9T 

CM2/SEC 

mobility  mu 
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6*6.9 
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EM  » 0/MU 
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EV 
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MESH  ■ 500.  CMAX  ■ 20.00  EV.  DC  • .0*0  EV. 


Fig.  4.12: 


Summary  of  plasma  parameters. 
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Fig.  4.13:  Tabular  electron  distribution 
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Fig.  4.18:  (Partial)  reproduction  of  tabular  summary  of 
reaction  sensitivities. 
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SUMMARY  OF  ELECTRICAL  and  optical  parameters 
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.08 

.03 

-.00 

-.00 

20.00 
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57.09 

.39 

6.00 

.16 

.00 

.00 

30.00 

5.69 

82.19 

2.38 

1.69 

.35 

.61 

.50 

00.00 

7.07 

102.29 

3.36 

6.75 

1.60 

6.16 

6.00 

50.00 

8.20 

119.09 

3.31 

5.52 

1.83 

9.7* 

8.06 

60.00 

8.70 

126.63 

3.53 
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2.65 
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60.00 
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90.00 

9.78 

1*2. S7 
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160.00 
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7.00 
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6.95 
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136.64 

7.61 
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220.00 

8.04 

129.71 

9.07 

23.22 

21.06 

9.69 

6.42 

230.00 

8.76 

128.32 

9.61 

26.09 

22.67 

9.66 

6.38 

260.00 

6.62 

126.76 

9.75 

26.93 

26.31 

9.57 

6.34 

250.00 

8.50 

125.02 

10.10 

2S.73 

26.00 

9.49 

6.29 

260.00 

6.31 

122.28 

10.69 

26.67 

27.76 

9.38 

6.25 

270.00 

8.13 

119.80 

10.90 

27.16 

29.62 

9.23 

6.18 

260.00 

7.97 

117.5* 

11.33 

27.92 

31.63 

9.09 

6.09 

290.00 

7.83 

115.51 

11.73 

28.72 

33.68 

8.96 

6.01 
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Fig.  4.21:  Tabular  summary  of  miscellaneous  electri- 
cal and  optical  parameters  as  a function  of 
time  (voltages,  current,  power  densities). 
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Fig.  4.22: 


Instantaneous  net  gain  coefficient,  rela- 
tive to  threshold,  as  a function  of  time. 
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Electrical  and  optical  power  densities. 
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4. 3 Error  Diagnostic  Assistance 

Although  it  is  not  possible  to  anticipate  all  possible  forms  of  incor- 
rect usage,  program  LASER  has  been  designed  to  provide  limited  di- 
agnostic assistance  upon  recognition  of  certain  (potential)  error  con- 
ditions. The  search  for  errors  occurs  in  processing  an  input  reaction 
queue  for  subroutine  synthesis,  and  at  the  time  of  data  input  for  execu- 
tion of  the  analysis. 

The  DATA  BLOCK  1 shown  in  Fig.  4.36  violates  several  of  the  conven- 
tions described  in  Sec.  3.2.  Assume  that  program  LASER  has  been 
dimensioned  NMAX  = 30,  KMAX  = 16,  and  NKMAX  = 8 (for  purposes 
of  illustration),  and  that  DATA  BLOCK  1 is  submitted  for  a program 
synthesis.  Even  if  the  reaction  syntax  contained  no  intrinsic  defects, 
its  species  content  and  total  length  would  exceed  the  dimension  capaci- 
ty specified  for  LASER. 

The  reaction  summary  generated  in  an  attempt  to  process  the  DATA 
BLOCK  1 of  Fig.  4.36  contains  numerous  informative  and  diagnostic 
comments,  as  shown  in  Fig.  4.37.  For  defective  reactions,  there  is 
no  FORTRAN  translation:  the  reaction  is  ignored  (and  unnumbered  in 
Fig.  4.37),  and  processing  continues.  The  data  file  generated  (unit  4) 
contains,  in  addition  to  the  rate  constant  data,  a logical  parameter 
(written  in  the  second  record)  whose  value  is  .TRUE,  or  .FALSE,  de- 
pending upon  whether  or  not  errors  were  encountered.  A later  at- 
tempt to  execute  an  incomplete  analysis  from  which  reactions  were 
rejected  can  be  prohibited  by  specifying  FATAL  (1)  = .TRUE. 

Return  now  to  the  KrF  program,  synthesized  from  the  reaction  scheme 
of  Fig.  4.1  to  produce  the  subroutines  illustrated  in  Fig.  4.2-  4.4.  For 
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the  purpose  of  illustrating  execution  error  diagnostics,  suppose  that 
the  original  kinetic  scheme  of  Fig.  4.1  had  contained  some  (unimpor- 
tant) error,  such  as  the  rejection  of  a duplicate  reaction.  (Thus,  for 
all  practical  purposes,  files  KRFBCD  and  KRFLASER  are  complete- 
ly satisfactory,  although  KRFDATA  will  contain  a record  of  the  fact 
that  a defect  had  been  encountered  during  synthesis. ) Suppose,  fur- 
thermore, that  the  dimension  declarators  in  LASER  have  been  changed 
since  the  original  creation  of  the  KrF  files,  and  that  execution  of  the 
analysis  is  attempted  with  the  input  DATA  BLOCK  4 shown  in  Fig.  4.38. 
Note  that  no  specification  of  the  vector  FATAL  (I)  appears  in  the  data 
group  $CONTROL  ...  $ to  override  the  default  initialization  .FALSE. 
Note  further  that  Fig.  4.38  differs  from  Fig.  4.  7 in  that  now,  there 
are  rate  constant  modifications  specified  by  $RATES  ...  $ entry,  and 
a new  species  "XXX"  has  been  initialized.  (Most  of  the  other  para- 
meters are  numerically  different,  but  that  is  not  relevant  for  the  point 
to  be  illustrated  here. ) Fig.  4.39  is  the  tabular  summary  of  the  reac- 
tion scheme  and  its  modified  rate  constants;  note  that  some  of  the 
attempted  modifications  were  illegal  and  were  rejected.  Fig.  4.40  is 
a summary  of  error  diagnostics  encountered  in  the  attempt  to  execute 
the  analysis. 
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Fig.  4.36:  Defective  DATA  BLOCK  1 (cont'd). 
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Fig.  4.39:  Tabular  summary  of  modified  rate  constants  for  the  kinetic  reaction  scheme. 
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Fig.  4.40:  Execution  error  diagnostics  illustrated. 
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